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FOREWORD 


This training guide explains the fundamentals and foundations of the long-range ballistic 
missile (LRBM) weapon systems launched from moving bases such as submarine or aircraft. 
Since the land-based missiles are essentially no different from the missiles launched from the 
moving bases except that the bases are fixed on the earth (without rotational and linear motions), 
the theory and analyses explained in this training guide would be equally applicable to all three 
classes of missiles—land-based, sea/undersea-based, or airborne missiles, which use inertial 
navigation system. 


It is written mainly for the training of new technical members using only introductory 
physics and calculus. Nonetheless, it manages to maintain the robustness of theory and the rigor 
of mathematical derivation. For this reason, it is recommended to veteran technical members as 
well, who might find a segment or two that would shed a new light on old subjects. 
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corrections; to Richard Hearn for typing unwieldly symbols; to W. H. Horton for his supervisory 
encouragement along with useful comments; and finally to author’s teacher and advisor, Professor 
Walter Hollister of MIT Department of Aeronautics and Astronautics for allowing him to use some 
figures on gyros from his book (Reference 1). 


This report was reviewed by William H. Horton, Head, Systems and Environment Branch 
of the SLBM Research and Analysis Division. 
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DAVID B. COLBY, Head 
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1.0 INTRODUCTION 


This publication is a training guide explaining fundamentals and foundations of the long- 
range ballistic missiles (LRBM) from the viewpoint of the navigation and guidance. The missiles 
may be launched from either fixed-base on land or from moving-base such as aircraft or 
submarine. The articles are written at such a level that anyone with freshman physics and calculus 
would not have much difficulty in understanding them. Readers of this publication are cautioned 
that the analysis of the actual, real-world systems are much more complicated than those described 
here. As the title implies, these articles are merely "fundamentals and foundations." 


Section 2.0 lists the symbols and notations used in this report. It includes an operator P ( -) 


for < (-). This is handy because the differentiation with respect to time in the earth-fixed E-frame 


may simply be expressed as P,. (-) . Section 3.0 describes the frames of reference relevant to this 


report. Section 4.0 explains implementation of Newton's second law of motion for the purpose of 
inertial navigation. Since all LRBMs use the center of the earth as the origin of the inertial 
reference frame, Section 5.0 deals with the errors incurred in the earth-centered inertial frame. 
Section 6.0 explains that the effect of angular rotations are dependent on the sequence in which 
they are executed, unlike vector additions. Theory of inertial navigation involves rotation of one 
reference frame relative to another. The theorem of Coriolis, which is discussed in Section 7.0, is 
essential to understand the relationship between the two rotating frames relative to each other. The 
matrix form of the theorem of Coriolis is discussed in Section 8.0. Section 9.0 discusses the 
rotational form of Newton's second law of motion, which is essential for the understandin g of the 
operation of gyros, gyro-accelerometers, and stable platforms of guidance system of LRBMs. 
Section 10.0 describes the prototype linear accelerometer, which serves as a tutorial for 
understanding the pulse integrated pendulous accelerometers (PIPA), and pendulous integrating 
gyro accelerometers (PIGAs), is described in Section 11.0. Section 12.0 discusses the rotational 
motion of the missile body based on some simple assumptions, which allows us to avoid 
complicated mathematics. 


Based on theoretical and mathematical tools developed in previous sections, the theory of 
practical gyros are explained in Section 13.0, the single degree-of-freedom (SDOF) gyro in Section 
14.0, the stable platform of inertial measurement unit in Section 15.0, and the PIGA in Section 
16.0. 


Moving from components to the system level, the derivation and implementation of the 
equation of motion for the inertial navigation is described in Section 17.0. The determination of 
the initial position and initial velocity to be used in the integrations by the guidance computer of the 
Newton's second law of motion is described in Section 18.0. Moving to the rocket steering 
scheme during powered stages, the so-called cross product steering is described in Section 19.0. 


Toward the end of the powered flight, the guidance system of the missile tries to correct its 
eITOF, in a Statistical sense, by viewing a predesignated star. The weighing matrix, also called the 
W matrix, which is used in the statistical computation for this purpose, is derived in Section 20.0. 
The topic of Schuler tuning is discussed concisely in Section 21.0. Schuler tuning is conceptually 
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important because it deals with the prevention of platform misorientation caused by the acceleration 
in the earth's gravitational field. Finally, the correlated velocity, which is the required velocity of 
the reentry vehicle to hit the target on the rotating earth by free falling with specified time-of-flight 
(TOF), is discussed in Section 22.0. | 
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2.0 LIST OF SYMBOLS AND NOTATIONS 


Rpg Displacement vector from point P to point Q. Similarly Vpg is the velocity vector of Q 
relative to P, etc. 


W,s Angular velocity of Frame B relative to Frame A. 


v“ Any vector, (e.g., R or W, etc.) with components expressed in Frame A. 
d 
P ee) = a= f + 
(=—0) 


=()= Ja 





P=) 


P, (-)=~(-) Time derivative with respect to (WRT ) Frame A, (ie., the increment observed in 


Frame A). 
Thus, 


A 


Xx 
PAR=P(R*)=<(R)=< y 
Z 


A 
2 
P?R = P,(P,R)=P(PR*)=5 


X 
dt?|” 
Zz 


P,P,R = P, (PR) 


The advantages of these notations will become obvious in the second-half of this report. 
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Let 


D D 
WwW 


x 


R? = and W,y =| w, | ; 


N < w 


WwW 


Zz 


D 


then the matrix equivalent of the vector W,,, 


denoted by W,, is given by 


0 -w, Ww, 
x 


DK_| w 0 -w 
Wis : 


Ww. OW 0 


so that [W>s [R= W>, xR? or 


O -w, wy \x i j —k W,Z— W.Y 
wz O -wWy ly|@iwy wy wz) | wW.x-w,z 
“Wy Wx 0 fz kX y Zh W,Y — WyX 


C. coordinate transformation matrix from Frame A to Frame B. 
Note: 
Ch3Cr =Cr 


CC, =C, =I identity matrix. 
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3.0 FRAMES OF REFERENCE 


3.1 INERTIAL (D-FRAME 


The frame in which Newton’s laws of motion (specifically the second law) is valid is called 
inertial (I)-frame by definition. In all, LRBM, an earth-centered nonrotating frame in the inertial 
space, is accepted as a (quasi) I-frame. Obviously, since the center of the earth accelerates around 
the sun in an elliptic orbit, the earth-centered I-frame (denoted by Ig) is not an I-frame in a strict 
sense. However, the error incurred in this practice is negligible (this is shown in detail under a 
separate heading in the sequel). 


In the standard I,-frame, the z-axis is collinear with the earth’s polar (spin) axis. The x- 
axis and y-axis are on the earth’s equatorial plane (but not rotating with the earth) in such a way as 
to form an orthogonal triad by the right-hand rule. This frame is sometimes referred to as the 
“polar earth-centered inertial frame." 


3.2 EARTH (E)-CENTERED, EARTH-FIXED FRAME 


This frame, being earth fixed, rotates with the earth relative to the inertial space. The x, y, 
and z axes are defined similar to the standard I, frame, except x and y axes are fixed to the earth 


and therefore rotate with the earth. 


3.3 NAVIGATIONAL (N), LOCAL-LEVEL FRAME 


The origin of this frame is located at the center of the onboard inertial navigation system 
(INS) of the ship, submarine, or aircraft cruising on or near the surface of the earth. The x-axis 
points to the north and y-axis to the east, both on the local-level plane, which is tangent to the 
reference ellipsoid. The z-axis points downward perpendicular to the tangential plane. The z-axis 
may or may not coincide with the local vertical (plumb bob) lim because the earth’s gravity field is 
not uniform. This causes the deflection of the vertical. 


3.4 GUIDANCE (G)-SYSTEM INERTIAL MEASUREMENT UNIT (IMU) FRAME 


The origin of the G-frame is located at the center of the missile guidance system’s IMU, 
with its axes parallel to the earth-centered I-frame. Obviously, since the missile is subject to 
acceleration, the G-frame cannot be an I-frame. However, since its axes are parallel to the axes of 


the I-frame, the coordinate transformation matrix from the G-frame to the I-frame denoted by Cc 


must be an identity matrix or C= I, where I denotes the identity matrix. 
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The onboard accelerometer measures the acceleration of the instrument relative to the origin 
of the I-frame, not relative to the origin of the G-frame to which the instrument has a fixed distance 
and therefore the relative acceleration should be zero. However, since the axes of G-frame is 
parallel to the axes of I-frame, the accelerometer output (sensed acceleration and the velocity-gain 
caused by it during the computation interval) should have the same values whether they are 
expressed in I-frame or G-frame. That is, 


AV' = CLAVS = IAV® = AV 
where 


AVS = velocity gain (AV) determined by the guidance system with components expressed in 
G-frame, and | 


AV! = AV with components expressed in I-frame. 





Thus 


AR’ =AV'At = AV°At 


1 


where 
AR! = position displacement with components expressed in the I-frame 


We see that as far as the determination of the missile’s incremental position and velocity are 
concerned the G-frame behaves as if it were an I-frame, although in reality it is not. To determine 
the current position or velocity relative to the I-frame, the incremental position or velocity obtained 
in the guidance frame is added by the guidance computer to the previous position or velocity, 
referenced to the origin of the I-frame, which is conventionally the center of the earth in all LRBM 
<r For this reason, some incorrectly call the G-frame “inertial,” which should be 
avoided. 
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4.0 NEWTON’S SECOND LAW 
(FOR THE IMPLEMENTATION OF INS) 


According to Newton’s second law 
mP’R, = LF | (4-1) 


Since the sum of externally applied forces XF has to be the sum of the gravitational force 
F, and the nongravitational force Fyg, we may write (1) as follows: 


2 _K Fro _ 
PrRyp = ae g+f (4-2) 


where g 1s the gravitational force per unit mass and f is the nongravitational force per unit mass. 


. By convention, f is called simply “specific force,” meaning eae (per unit mass) 
nongravitational force. 


2 
The reason we separate f=—— = from g=—— 1s by necessity, not by choice, because the 


on-board inertial instrument lisesiliiitsiiins can measure only f, not g. This is explained in the 
section of the accelerometer later. The g is determined by the onboard computer, based on the 
mathematical model (inverse square law) as a function of the position Ry from the center of the 
earth, which is the origin of the I-frame by convention. 


Thus, in the INS, Newton’s second law (4-2) is implemented as shown in Figure 4-1. 





me 
iN 
N < x< 


position 







SERVO-CONTROLLED 
INITIALLY NON- 
ROTATING PLATFORM 


FIGURE 4-1. NEWTON’S SECOND LAW IMPLEMENTATION OF INS 


rs 
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5.0 ERRORS IN THE EARTH-CENTERED INERTIAL FRAME (Ig FRAME) 


For the purpose of estimating the magnitude of the errors incurred by treating the earth- 
centered inertially nonrotating frame (I; frame) as the I-frame, consider a system consisting of the 
sun (S), the earth (E), the moon (M), and a missile (P) near the earth with the center of the sun 
considered as the origin of an I-frame as shown in Figure 5-1 (Figure not to scale; the rest of the 
universe is ignored). 


E 


FIGURE 5-1. THE SUN, MOON, AND EARTH SYSTEMS 


The forces acting on missile P with mass m are the gravitational forces by the sun, the 
earth, the moon, and the nongravitational forces of thrust and aerodynamic force (if in the 
atmosphere). 


By applying Newton’s second law along with the law of gravitation to the missile at P. 
_GM,m 


GM,m GM,,m 








RARE, eee Fo = mi Be (5-1) 

m = _ mass of the missile at P 

M, = mass of the sunatS 

Me = mass of the earth atE 

Mm = mass of the moon atM 

Usp = unit vector along SP 

Uep = unit vector along EP 

Ump = unit vector along MP 

G = gravitational constant 
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Since Resp = Rs + Rep 
mPRgp ae mP; Re + mP; Rep _ (5-2) 
Dividing by (1), by m and using (2): 


GM Ugg — Me Ugg ~My +B = PPRgp + PPR 
R¢p Rip Rep (5-3) 











Now, applying Newton’s second law to the earth: 








_GM,M, GM,,M, 
M,PrR, = R? ——— U,, = R?, —_— —— Use 
SB (5-4) 
or 
GM GM 
PPR, =-— =U MU ve - 
1 *SsE R2, Use — R?, UME (5-5) 
Substituting (5) into (3) for PIR = 
P’Rep —_ ie Usp 4 Fro Ags + Agy 


Rip )StiéiaM (5-6) 


where 


1 1 
Ag. = —GM,| ——-U;p -—- U 
. {ge RR «| (5-7) 


is the difference between gravitational forces per unit mass by the sun on the missile and by the sun 
on the earth, and 


1 1 
Ag,, = -GM Ui .== =U 
. dg oy Rare | (5-8) 


is the difference between gravitional forces per unit mass by the moon on the missile and that by 
the moon on the earth. 


For a missile about 1,000 km above the surface of the earth: 
Ag, 
—— is less than one part per 10 millions in magnitude. 


E 
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Agy 





is less than two parts per 10 millions in magnitude, 
E 


GM 
R? = is the magnitude of the earth’s gravitational force on the missile per unit mass. 
a | 





where g.= 


These values are at least one order of magnitude smaller than the accuracy of the best state- 
of-the-art accelerometer even in the laboratory environment. 


Although we have not considered the effect of the rest of the universe on the missile, it 
turns out that its magnitude is much smaller than the effects from the sun and moon. 


For these reasons, (6) is approximated by 





a Up + 


Se 
, a 


PrRpp =— 
ia a m (9) 


with negligible error represented by (7) and (8). 


| =F 
Equation (9) is the exact form of Newton’s second law SF = ma or a= = as if the 


origin of the I-frame were located at the center of the earth. 


Qualitatively speaking, the magnitude of Ag, as a proportion in Pr R,»p 1s neglible because 


the sun pulls both the earth and missile together so that the relative position of the missile with 
respect to the earth is relatively unaffected. Similar reasoning applies to the case of Ag,, being 


negligible. 
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6.0 ANGULAR ROTATIONS 


In the addition of two vectors V; and V2, V; + V2 = V2 + V; as shown in Figure 6-1, 
(i.e., Vy followed by V2 has the same result as V2 followed by V}). 


V1 


V2 V2 


V1 


FIGURE 6-1. VECTOR ADDITION 
In angular rotations, a rotation about the x-axes followed by another rotation about the 
displaced y-axes does not result in the same orientation of the frame as a rotation about the y-axes 
followed by another rotation about the displaced x-axes. In this sense, angular rotations do not 
commute, (1.e., do not act like vectors). 


To see this by way of two simple demonstrations, consider an orthogonal (Frame-0) with 
Xo, Yo, and Zp axes, as shown in Figure 6-2. 


Frame 0 


Zo 


Yo Yo 
XO 


Xo 


FIGURE 6-2, O-FRAME 


6-1 
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In the first demonstration, rotate the above frame about XQ - axes by 90° (using the right- 
hand rule). This results in Frame 1 (with X;, Y;, and Z; axes) as shown in Figure 6-3. 


Frame 1 


Zi Zi 
X1 = Xo 


X1= Xo 


FIGURE 6-3. FRAME 1 


Next, rotate Frame-1 about Y; - axes by 90°—esulting in Frame-2 (with X2, Y2, Zz axes) 
as shown in Figure 6-4. 


Frame 2 


Y2= Yi Y2=Yi 


FIGURE 6-4. FRAME 2 


Now in the second demonstration, returning to Figure 6-2, rotate Frame 0 about Yo-axes 
(instead of Xo axes) by 90°—resulting in Frame 1° (with X;°, Y;°, and Z;° axes) as shown in 
Figure 6-5. 


6-2 
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Frame 1° 


Z1 
Y*°1 = Yo Y‘1= Yo 
Z1 
> x) 


FIGURE 6-5. FRAME 1° 


Next, rotate Frame 1° about X,°-axes by 90°—resulting in Frame 2° (with X2°, Y>’, and 


Z2° axes) as shown in Figure 6-6. 





Frame 2” 


Z‘2 Z‘2 Y°2 
K 
Y2 


X2=X1 X2= X11 
FIGURE 6-6. FRAME 2’ 
Obviously, Frame 2° (Figure 6-6) has a different orientation from Frame 2 (Figure 6-4). It 
means that, in describing the effects of angular rotations, it is necessary to specify the order in 


which a sequence of rotations is executed. 


Consider Frames A, B, C, and a vector R, which has components RA = ( Xa, Ya, Za )T in 
Frame A, RB = ( Xp, Ypg, Zp )T in Frame B, and RC = ( Xc, Yc, Zc )T in Frame C. 


Then 
R* = CR“ (6-1) 


and 
R© = CER® (6-2) 


Substituting (6-1) into (6-2) for RB 
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R°=C,C,R* (6-3) 
Since 

R°=C,R*. (6-4) 
From (6-3) and (6-4), we conclude | 

CEC, = Cy. (6-5) 


Referring to (6-5) and based on our discussions, C;C;, #C, C, 


As a matter of fact, C. c- has neither physical nor geometrical meaning. In our notation, 
the second rotation C. is written to the left of the first rotation Cc, so that the superscript B of Cc, 
matches the subscript B of C5, so that C,C, becomes C; as in (6-5). 


We can extend this process, for instance, to: 
CCC C.-C. (6-6) 
By writing as in (6-5) or (6-6), we can check the consistency of multiplications of the 


coordinate transformation matrices. If superscripts and subscripts do not match in a way we have 
specified, it indicates errors somewhere. 
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7.0 THEOREM OF CORIOLIS 


Consider two references, Frames A and B, with common origin at 0. Frame B rotates with 
respect to Frame A with angular velocity Wap. A vector Rop = R may be expressed in Frame A 


R* =IX+JY+KZ (7-1) 


and in Frame B 


R® =ix+jyt+kz. (7-2) 
Now, using (7-1): 


dX _dY _dZ 
=f 4 kK 
dt “dt dt (7-3) 


because the unit vectors I, J, and K are fixed in Frame A, and do not vary with time in Frame A. 


Using (7-2): 


dx .dy .dz di dj dk 
jp ek ey 
dt ‘dt dt 7 dt Yy dt * dt (7-4) 


Note: In this case, the unit vectors i, j, and k, which are fixed in Frame B, rotate relative to 
Frame A, and therefore, are time-varying as viewed from Frame A. 


Now, consider di in (7-4). (Referring to Figure 7-1). 
dt 


Since vector i is a unit vector, it cannot change its magnitude. However, it can change its 
direction relative to Frame A because 1 is fixed in Frame B, which rotates relative to Frame A. For 
an infinitesimal angular displacement, we can visualize the tip of the i-vector moving on the plane 
parallel to the j - k plane. Therefore, we may express Ai caused by the angular displacement of i in 
Frame A in terms of its displacements in the j and k directions. 


Referring to the Figure 7-1 below, we see 
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W, 





FIGURE 7-1. INCREMENTAL ROTATION OF FRAME 


that W, causes Ai = jW,At in the j direction during At and W, causes Ai = -kW,At during At in the 
—k direction. Summing the two components, we have 


Ai = jW_At - kW,At (7-5) 
dividing by At and taking the limit, we have 
di 


— = jw. -kWw 
Loo (7-6) 


It turns out that the right-hand side of (7-6) is also equal to WB xi as shown below: 


= jW,-kWy (7-7) 


7-2 
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It follows: 
oe wx i 
dt (7-8) 
Similarly, we can show that: 
ae Ww xj (7-9) 
—=w xk. 
dt (7-10) 
Now, referring to the right-hand side of (7-4), and using (7-8), (7-9), and (7-10): 
di dj dk 
Xx—+ y—+zZ— 
dt “dt dt 
=X (w? xi)+y(w? xj)+z(w? xk) 
=w” x (ix+jyt+kz) 
=w’ xix+w” xjy+w, xkz (since x, y, and z are scalars) 
=w'xR® 
= Was xR? (since W = Wag by definition in our case) (11) 
Now, 
p,R =P (ixt+jy+kz)” 
_;% oy 
dt : dt dt (12) 
because 1, j, k are fixed in Frame B. 
So, substituting (7-12) and (7-11) into (7-4): 
P,R® = P,R® + Wry XR?. (13) 


7-3 
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In (13), each term may be expressed, after appropriate transformations, in either Frame A 
or Frame B, (see next section). So, dropping the superscript B from (7-13), 


P,R = PaR+WapxR (7-14) 
which is called the Equation of Coriolis. 


What the Equation of Coriolis implies is that the differentiation (WRT time) of a vector in 
one frame is not equal to the differentiation of the same vector in another frame, which is rotating 
with respect to the first frame. And to obtain the value of PR of (7-3) in terms of PgR of (7-12), 
we have to add a correction term W ap x R (which incorporates the effect of rotation of Frame B 
relative to Frame A) to PgR as shown in (7-14). 


For some, the Equation of Coriolis is considered the second most fundamental equation 
only next to the Newton’s second law in the analysis of navigation and guidance. Readers who are 
interested in the geometrical approach in the derivation of (7-14) may find it in other text books 
such as "Mechanics" by Symon (published by Addison Wesley). Some may find the geometrical 
approach simpler and easier to follow, while others may not. For this reason, an analytic approach 
is presented here to assist the comprehension in view of the importance of the theorem. 


7-4 
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8.0 MATRIX FORM OF THE THEOREM OF CORIOLIS 


The law of Coriolis in vector form (which we have derived previously) between Frames A 
and B follows: 


P,R = P3R + Wag XR. (8-1) 
Coordinating (expressing components) in Frame B: 
(P,R) =(P,R) + W2, xR®. (8-2) 
Now, referring to (8-2): 


(P,R) =Ci(P,R)" 
= C®(PR“) because (P,R)* = (PR*)" = PR“ 
= C,[P(C3R’)| 
=C®| (Pct )R®+C4 (PR®)] 
= C(PC*)R® + ChCPR® 
= C(PC*)R® + PR® (8-3) 


. BRA B 
since C, Cp =C =I 
Returning to (8-2), recall the following from Section 1: 
Was XR® «> [Way |R” (8-4) 


where for Wis = (WW, W, \" 


8-1 
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0 -W, W, 
Wk=| W, 0 -W, 
-w, Ww. 0 


which is the matrix representation of the vector cross—product. 


Substituting (8-3) and (8-4) into (8-2) 


C?(PC4)R® +PR® = PR" +[W2]R® om, 


noting that (P,R) =(PR®) =PR®. 
It follows: 


Ci(PC3)R” =[Was JR”, | (8-6) 


Since (8-6) has to be valid for all R®, we have 
C,P(C3)= Was (8-7) 
premultiplying both sides of (8-7) by 
(Ch) =(Cx) = CS 
because Cc, , being a coordinate transformation matrix (between two orthogonal frames), is an 


orthogonal matrix, we have 


CECyP(Cz) =CiWas (8-8) 
since C; Ch = Cc; =| 
PC} =CiW.. (8-9) 


Equation (8-9) is a matrix form of the Equation of Coriolis. 
Since 


,_d Cé(n+1)—C4(n) 
PC) = —C, = + 
» dt At (8-10) 


8-2 
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where C5 (n) denotes the value of C} at nT. 


We have from (8-9): 
Cp (n +1) =C3(n)+C3(n)Wy (n)At - (8-11) 


Equation (8-10) suggests that if we can measure Wis by gyros fixed to the body (such as 
the rocket structure), then we can update the orientation of the body-fixed frame relative to a 
reference frame (such as the I-frame) by updating the Cy matrix. This idea is used in the so-called 


strap-down INS where the gyros and accelerometers are fixed to the body instead of on the inertial 
platform, which is isolated from the body motion by a servo system with gimbals. It may also be 
used in error analysis of the IMU platform drift caused by the gyro-drifts rates. 


8-3 
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9.0 ROTATIONAL FORM OF NEWTON’S SECOND LAW 


The rotational form of Newton’s second law is the basic starting equation essential to 
describing the operation of gyro, gyro-accelerometer, and base—motion isolation systems. One 
such system is the stable platform of the IMU of a LRBM, which is controlled by a servo- 
mechanism using gyros and gimbals. 


The equation also leads to the so-called Euler equation, describing the dynamic behavior of 
the rigid body motion, which is used in the SDOF simulation model of the rocket body in flight. 


Consider a point mass m located at point k. By definition the moment of momentum or 
“angular momentum” Hy, of m about the point I (origin of the I-frame) is 


H,,=R,xmP,R, . (9-1) 


The moment of force or “torque” My, of external force on m at k about the point I is by definition 
the following: 


M, =R, xK (9-2) 


Differentiating (9-1) WRT time and noting that the cross-product of two identical vectors is 
Zero: 


PH, =P,(R, x mPR,) 


=PR, xmPR, +R, xmPR, 
=R, xmP’R, 
=R, XF, 
It follows in view of (9-2): 
PH, = My, (9-3) 


Summing up (9-3) for all particles (all k’s), and defining Hy and My by 


H, = D,(R, <m,PR, ) (9-4) 


M, = ¥.(R, xF,) (9-5) 


and with the understanding that F, refers to the external forces only, we get 
PH, =M,. (9-6) 
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This is the rotational form of Newton’s second law. The law states that the time rate of 
change (WRT the I-frame) of the angular momentum about an inertial point is equal to the torque 
about the same inertial point. | 


By denoting the location of the center of mass by C, Ry. may be expressed by a simple 
vector addition: 


R,=R,.+Re (9-7) 
Substituting (9-7) into (9-4) and (9-5) and then the results into (9-6), after some algebra 


(see for instance Reference (9-1)), we get the rotational form of Newton’s second law about the 
center of mass (instead of the inertial point): 


PH, =M, | (9-8) 
where 
H.=2, (Ro, xm,PRy, | (9-9) 


= the angular momentum of the system of particles about the center of mass 
M, = >. Ra, x F, (9-10) 


= the external torque about the center of mass. 


Equation (9-8) states that the applied torque about the center of mass is equal to the time 
rate of change (WRT the I-frame) of the angular momentum about the center of mass. 


The rotational form of Newton’s second law is valid for any system of particles whose 
internal forces are central, whether or not the body is rigid. 


Pa) Swoweerry cee attstiinnn eee —nncctiecemenl, a sthtmentatmemmten ee) eae cee een eee. et eee ll sears mien = =©6aimat 
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10. PROTOTYPE LINEAR ACCELEROMETER* 


Consider a box (instrument) with a spring and a test mass m. The box is rigidly attached to 
the rocket structure as shown in Figure 10-1. The rocket is moving along a straight line under a 
constant acceleration. 


Under this steady-state condition, the spring extends by a constant displacement x from P 
to Q along the unit vector U, because of the inertial reaction force on m, which is the reaction to the 
applied force exerted on P. 


Direction of 
the total 
/ acceleration 


I-frame origin 





FIGURE 10-1. PROTOTYPE LINEAR ACCELEROMETER 
The forces acting on m at point Q are the gravitational forces Fgg (the gravitational force Fg 
at point Q), and the force -kx pulling the mass toward the point P by the restoring force of the 
spring. 


10-1 
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Applying Newton’s second law to m at Q with components expressed along the missile 
center line, by taking dot product with u: 


mP?Rg -U = (Fog — kx): 0 | (10-1) 


Since by vector additions: 
RY =R,+R, (10-2) 
we get by substituting (10-2) into (10-1): 


(mP’R, + mP/Rjg)-u = (Fog — kx): u (10-3) 


Under steady-state conditions Rpg is constant, which makes PjRpg and, therefore, P? Royo 
equal to zero. It follows from (10-3): 


F 
PR,-u=(—2-Ex)-y 


m m- (10-4) 


F 
Noting that a is the gravitational acceleration g at Q (denoted by go), we have from 
(10-4): 


2 k 
PPRy-u=(gq~ =x) 0. (10-5) 


Applying Newton’s second law to the whole box considered as a rigid body (instrument with the 
pivot point at P ), noting that the acceleration at P has to be caused by either gravitational specific 
force gp at P or nongravitational specific force fyGp at P, and taking the component along the u 
direction: 

(PrR,) ene (g, + fice) "u (10-6) 


Equating (10-5) to (10-6) and rearranging 


f w= xt(g —g )-u 
NGP m Q P (10-7) 


10-2 
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Now both gp and gg are gravitational accelerations (gravitational force per unit mass) several 


thousand miles away from the center of the earth at points P and Q that are less than only a few 
centimeters apart. Thus, gq and gp are essentially equal and go - gp = 0. 


It follows from (10-7) 


f, we 
Nop (10-8) 


Equation (10-8) shows that the mass-spring box (an idealized accelerometer) measures the 
externally applied nongravitational force per unit mass (called nonfield “specific force”), which is 
equal to the force (per unit mass) that the test mass m exerts on its support at P through the spring. 
That is, the spring displacement x is the measure of nongravitational force per unit mass (specific 
force) along the direction u at point P. Since x is preportional to -fygp - u, we should say, more 
exactly, that the accelerometer measures the negative of the nonfield specific force along (u). 


It is important to note that the accelerometer (a misnomer in a strict sense) measures 
nongravitational specific force fg, not fyg + fg , Newton’s second law states that 


PP Rp = of = fg + fig - 


Thus, in the implementation of the navigation system that relies on the onboard 
accelerometers, we have to always add the gravitational acceleration (computed generally based on the 
inverse-square gravity model) to obtain the total acceleration as shown below in Figure 10-2. 


—____ : 
accelerometer total acceleration 


f, 


G - computer 
FIGURE 10-2. COMPUTATION OF TOTAL ACCELERATION 


For some, the fact that the accelerometer onboard a rocket in motion under a gravitational 
field measures only the nongravitational specific force may seem counter-intuitive. 


| Referring to the first figure of this section, both points P and Q located more than 6,000 km 
away from the center of the earth are only less than a few centimeters apart. This means that the 
gravitational force per unit mass at these two points is essentially equal. 
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Thus, we can see that the effect of gravitational force on the displacement of the spring is 
essentially zero. In another words, the accelerometer is incapable of sensing the gravitational force 
and therefore cannot measure it. Now, let us consider the combined effect of thrust and aero- 
dynamic forces along the direction u. This will push the rocket structure forward along with the 
pivot-point P, which is fixed to the structure. By Newton’s first law, the mass M located at Q will 
resist this change and try to remain behind, thus stretching the spring by inertial reaction. Now, 
we can see why the accelerometer senses and measures the effect of nongravitational force only 
excluding that of gravitational force. | 


10-4 
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11.0 PULSE INTEGRATED PENDULOUS ACCELEROMETER (PIPA)* 


Consider a physical pendulum of mass m with the moment of inertia J and with the center 
of mass (CM) located at distance L from the pivot point P. It is damped with the damping 
coefficient C. | 


The vehicle with the pendulum (pivot point P) fixed to it is under acceleration by the 
combined effect of both gravitational and nongravitational forces. The nongravitational portion of 
the acceleration is acting in the direction labeled as IA in Figure 11-1. 

F, 
As explained in the case of linear accelerometer, the gravitational force per unit mass aa te 


or gravitational acceleration at pivot point is essentially identical to that at the CM of the pendulum. 


P 





Reference Line 


FIGURE 11-1. PHYSICAL PENDULUM 


Thus, the gravitational acceleration does not contribute to the swinging of the pendulum. 
However, the nongravitational specific force fg (thrust and aerodynamic force) is applied to the 
point p, which is fixed to the vehicle, but not to CM (unlike the gravitational force that is exerted on 
both p and CM). As p accelerates along IA, m, according to Newton’s first law, tries to remain 
Stationary in the inertial space, thus making angle 6 with the reference line. 


Applying the rotational form of Newton’s second law about a pivot axis (into the paper): 


Jo=>M (11-1) 


*This section is adapted from Chapter 1.0 (the introduction) of MIT Instrumentation Laboratory Report R-574 


(on PIPA) by Kee Soon Chun, March, 1967. 


11-1 
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The applied torques 1M consists of the following three parts: 


F 
¢Torque caused by the inertial reaction force —Fyg=—M— =~ mfg about IA axis. 


The negative sign is necessary because the inertial reaction force is in the opposite 
direction from the applied Fyg. This torque is from the diagram, Lcos@(-mfy,)=- 
mLcosOfnxg 


«Damping (dynamic friction) torque - C8 (with negative sign because it opposes the 
motion). 
*Rebalancing torque generated by the torque generator Mig, which drives the angular 


deviation 9 to zero. 


In actual operation, the angular deviation 6 of the pendulum from the reference line is 
limited to very small angles and constantly rebalanced by the rebalancing torque generated by the 
torque generator located at the pivot point (not shown in the diagram). 


Thus, from (11-1), summing up torques: 
J@= M,,— mLcos@f,,,—- C (11-2) 
since @ is limited to very small angles cos8 = 1. Thus 


J6+CO=M,-mLf,, (11-3) 
Integrating with respect to time for the duration of At: 


At At At At 


1 {6a+c | édt= fM_dt-mb | f,,ct 
0 0 0 0 (11-4) 


Because of the nature of the rebalancing operation, the angular displacement constantly 


reverses its direction. This causes the average value of @ (and average value of 6 as well) to 
approach zero for the time duration much greater than the period Ti, of each rebalancing torque 
pulse, which is indeed the case in actual application. 


Thus, for At >> ie 


[éat = fat =( (11-5) 


a OE 


, 
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It follows from (11-4) with (11-5): 


At At 
I, M,,dt = mL | fygdt . : (1 1-6) 


In the actual system, Mi, is proportional to the constant DC current I, through the torque 
generator (with the permanent magnet). Thus, 


Meg = Stg I | (1 1-7) 


where Sy, is the torque generator sensitivity. It follows, 
At At 
J, Mudt= | S,ldt=S,1At. (11-8) 


By making the integration duration At sufficiently small ( while keeping At >> T;, ) relative 
to the thrust and aero—dynamic force variation, we may treat f,, constant during At. Thus, 


At 
[ fugdt = fy At = AVyue ‘ink 


where AVng is the velocity increment caused by fg during At. Substituting (11-8) and (11-9) into 
(11-6): 


SiglAt — mMLAVxg 
Or 


S 
—_ _ 8 
AVyg =—£IAt. bial 


Thus, AVpg is proportional to I At (integration of current), or the amount of electric charge 
passed through the DC torque motor during At. If we use electric pulses 


with constant magnitude (instead of DC current), the summation (integration) of pulses is also 
equal to the amount of electric charges. Because AVng is the integration of fy, (which, in turn, 


proportional to the summation (integration) of current pulses) during At, this type of instrument is 
called pulse integrated. pendulous accelerometer (PIPA). 


11-3 


From (11-10) with At fixed 


' OF 
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S_At 

d(AV_, ) = —*—dlI 

(AV) mL 


AV, (using (10) again) 


om NG dl] 


I 


d(AV,,, ) 


AV xe 





(11-11) 


(11-12) 


Equation (11-12) shows that one part per million error in the DC current measurement 
would cause the corresponding one part per million error in the estimation of the velocity increment 
AVng. In actual implementation, the constant current source for the DC motor is achieved, in 
essence, by applying a constant voltage from a zener diode across a precision resister. But, both 
zener diodes and precision resistors are not only temperature sensitive but also deteriorate with 
time. For this reason, to achieve more accuracy, PIPA is replaced by PIGA in the newer long- 
range ballistic missile. PIGA is discussed in a later section. 


11-4 
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12.0 ROTATIONAL MOTION OF THE MISSILE BODY 
In the previous section on the rotational term of Newton’s second law (Section 9.0), the 
following equation has been derived: 
PiHc = Mc (12-1) 
where H, is the angular momentum about the center of mass, M, is the applied moment (torque) 


about the center of mass, and P;(-) is = (-) in the I-frame. 


Applying the theorem of Coriolis to (12-1) between the I-frame and a body fixed body (B)- 
frame of the missile and expressing the components in the B-frame: 


(P.H,) =P,H8 + We x He = ME (12-2) 


As shown in Section 9.0, the angular momentum H, is given by 


H.= dm,Rex x PReg jibes 


Applying the theorem of Coriolis to R,x between the I-frame and the B-frame, which 
rotates relative to the I-frame with the angular velocity Wyp, we have: 


PiRcx = PpRcx + Wy < Rex 

= Wipe X Rex 
since PgRcx =0 because Rcx is fixed in the B-frame. (Strictly speaking, rocket exhaust gases 
cause the center of mass of the rocket to shift, and Rex 1s not fixed in the B-frame. The 


derivation here assumes that the shift of the center of mass 1s negligible for the time span of this 
discussion. ) 


Substituting Wyg x Rcx for PyRcx in (12-3), we have 


H =) m,R,, x (W,, xR) (12-3b) 
K 


Using a vector identity: 


V,x(V,xV,)=(V,°V,)V,-(V,°V,)V 


3 


(V,-V,)V,-(V,-V,)V, for Vi=V3 
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we get from (12-3b): 


H.= Ym, (Rex Ra)Wp = Ym, (Rex Wa )Rex . | (12-4) 


At this point, we define the x, y, z axis of the B-frame with the origin at center of gravity (CG) of 
the missile as shown in Figure 12-1. 





FIGURE 12-1. B-FRAME FIXED TO THE MISSILE 


12-2 
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where the axes are defined as : 
X axis— Along the missile center line 
Y axis —Perpendicular to the center line at the CG 


Z axis —Perpendicular to both the x and y axis at the origin so as to form an orthogonal 
triad by the right-hand rule (X © Y = Z) 


Now consider the transversal cross-section of the missile on the y-z plane as shown below in 
Figure 12-2. 


FIGURE 12-2. CROSS-SECTION (y-z PLANE) OF THE MISSILE 


We see that for any point P}(y},Z}), there exist Po(y2,z2),within the missile such that 
YiZ, t y.Z, =ab+(—ab)=0. 


Extending this situation to cover the entire volume of the missile, 


=0 
DYn% (12-5) 


12-3 
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Next, consider a longitudinal cross-section on the x-y plane as shown below in Figure 12-3. 





FIGURE 12-3. CROSS-SECTION (x - y PLANE) OF THE MISSILE 


We see that for any point P3(x3, y3), there exist P4(x4, y4), within the missile such that 
X3V3 +X,y, =cd+(—cd)=0, 


Therefore, for the entire volume of the missile, 


X,Vx =O. 
2 Bue (12-6) 
Similarly on the xz plane: 
=0. 
| 2x7 (12-7) 
By denoting R¢, and W,, in (12-4) by 
12-4 
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Ky W, 
Ro. =| y, | and W, =| W 
CK K 1B y 
Zx W, (12-8) 
we have 
WwW, 
B pB B 24.2452 \|W 
ym, (RE, Re, |We=) my, (xe+y2+ 22 y (12-9) 
k Ww 


Zz 


Using (12-5), (12-6), and (12-7) in the second term of (12-4): 


Xy 
>> m, (Rex Wa JRex = ym, (x, Wx +y, Wy + Zx Wz) Yx 
K K Zy 
ym, x7 W, 
= ym y<W, y 
> m,z,W, (12-10) 
Substituting (12-9) and (12-10) into (12-4): 
ym, (yz + ZK )w, J, W, 
Hé =| >) m, (xz +25 )W, |=| JW, ~ (12-11) 
yim, (xz + Vx JW, JW, 
which may be written as : 
J, O O)/(W, 
He=|O J, O| |W, |=Jw (12-12) 
O O J) \W, 
where 
= >mx(yg +Zx) I 0 0 
y= > mg(xe +z) andJ=/0 J, 0}. (12-13) 
0 J 
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In (12-12) and (12-13), we note that all off-diagonal elements of the J-matrix, which is 
called the moment of inertia tensor, are all zeros, that is, the J-matrix is a diagonal matrix, thus 
greatly simplifying analysis and computations. The choice of body axes, which result in the 
diagonal J-matrix is called the "principal axes." What seems to be a natural choice of axes in our 
case happens to coincide with the principle axes. 


For analytical methods (which are not simple) for determining the principal axes, refer to 
Mechanics by K.R. Symon or Classical Mechanics by H. Goldstein. 


Two proven theorems for selecting the principal axes follow: 

¢ Any axis of symmetry of a body is a principal axis such as the x-axis in our case. 

¢ Any plane of symmetry of a body is perpendicular to a principal axis. In our case, the 
plane (containing the x-axis), which is perpendicular to the y-axis is perpendicular to the 
principal axis, thus making the y-axis the principal axis. Similarly for the z-axis. 


Now referring to (12-2), we have, using (12-11): 


d JW, | rw, 
P,H¢ = —| J,W, |= I,W, (12-14) 
LW, IW, 


since J, Jy, Jz are constants in view of (12-13) for constant fixed mass missiles. (Jx, Jy, Jz are 
not constants for variable mass missiles such as real-world missiles in flight.) 


And, 
i j k 
WexHo=|W, W, W, 
IW, J,W, JW, 
J,W,W, -J,W,W,) ((J,—-J,)W,W, 
=! J,W.W, -J,W,W, |=| (J, -J,)W,W, |. (12-15) 
J,W,W, -J,W,W,) ((J, -J,)W,W, 


Substituting (12-14) and (12-15) into (12-2), and expressing the components of M- by Mx, My, 
M, we have 


12-6 
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I,W, +(J,—J,)W,W, =M, (12-16a) 
J,W, +(J, —J,)W,W, =M, (12-16b) 
J,W, +(J, —J,)W,W, =M, . (12-166) 


Equation (12-16) is called Euler’s equations (of motion) for a rigid body (along the 
principal axes). : 


Equation (12-16) may be solved for W,(t), W(t), W.(t), (which are the angular velocities 
of the missile body relative to the inertial space) in terms of M,, My, M;. Mx, My, and M; are the 
externally applied torques generated by the reaction force to thrust (which is, by custom, simply 
called thrust) and the aerodynamic force about the CG of the missile (known as a function of time 
in simulation). 


In Section 17.0 on the matrix form of the theorem of Coriolis, we have derived 


d 


ae =C,Wa 


(12-17) 


(with the superscript A of (12-9) in that section replaced by I) where W;, is a 3X3 matrix given by 


O -wW, W, 
Ws =| W, O -W. 
-W, W, O (12-18) 


With W,(t), Wy(t), W,(t) from (12-18), then Equation (12-17) may be solved for Ci (t) 


Ci,(t) Cyp(t) C(t) 
C(t) =| C.,(t) C,(t) C.,(t) 
C5(t) C(t) C,,(t) (12-19) 


which is the coordinate transformation matrix from the B-frame to the I-frame, thus specifying the 
angular orientation of the missile in the inertial space. 
Newton’s second law applied to the CG of the missile and expressed component-wise in 
the inertial frame is: 
PPX, = fcx + fox (12-20a) 


P ry, = foy + fray (12-20b) 


12-7 
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PPZ, = fgz + faz (12-20c) 


2 
where Pe ( -) 7 (-) in the inertial space, f, is the gravitational force per unit mass, and fng is 
t 
the nongravitational specific (per unit mass) force which includes both thrust and aerodynamic 


forces. The fxg is most conveniently expressed in the B-frame and transformed to the inertial 
frame. 


The three equations of (12-20) completely specify the translational motion (position and 
velocity) of the missile’s CG along the inertial X;, Y], Zr axes, while another three equations 


(12-16) completely specify the angular orientation (which is the integration of the angular rates) of. 


the missile in the inertial space. For this reason, the mathematical model for simulation based on 
these six equations are called a six degree-of-freedom (6 DOF) model. 


If we are concerned only with the translational motion of the CG, ignoring the missile 
orientation, only the three equations of 20’s are used, and its model is called a 3 DOF model. 


Since the motion of CG of the missile is the same as the motion of the point mass in which 
the entire mass of the missile is concentrated at the CG, it is also called the point mass model. 


12-8 
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13.0 THE GYROS* 


The rotational form of Newton’s second law about the center of mass C is, as derived 
previously : 


PiH¢ = Mc. (13-1) 


GYRO ELEMENT GIMBAL 


Yge 
ae 


GYRO ELEMENT 
GIMBAL AXIS 





FIGURE 13-1. GYRO ELEMENT 
The assembly (the spinning rotor and its support gimbal) shown in Figure 13-1 is called 
“gyro element” (ge), with axis of the ge frame as shown in the figure. Replacing subscript C by ge 
(cm of ge) in (13-1) and, applying the law of Coriolis between I-frame and ge frame, 
PyHge = PpceHge + Wi-ge X Hee = Mge (13-2) 
where the subscript i-ge is the same as I-ge. 


The angular momentum of the gyro element Hg, consists of the spin part H, (constant) and 
nonspin part Hyg, that is : 


Hge = H, + Hps (13-3) 
where Hs is the angular momentum vector of the gyro element other than that of the spin rotor. 


Substituting (13-3) into (13-2): 


13-1 
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PyeHs + PgeHns + Wi-ge X (Hs + Hns) =Mge.  ~ (13-4) 


In (13-4), PseH; = 0 because H, is constant and fixed in ge frame, and therefore 


Wize x (H; + Has) = Wie X Hi; 


because 
H, >> Uys . 


It follows that 
PseHns + Wi-ze x H,; = Mge ° (13-5) 


‘When P,.Hns is negligible compared to Wj.ge X Hs, (13-5) reduces to 


Wize X Hy = Mge (13-6) 


Equation (13-6) is the approximate performance equation for the practical gyro. The 
geometric interpretation of (13-6) is depicted in Figure 13-2 with axes defined in Figure 13-3. 
That is, the spin angular momentum H, precesses (rotates) relative to the inertial space with the 
angular velocity of ge (relative to I-frame), Wj.ze, in an attempt to align itself with the applied 
torque vector Mge about the center of mass of ge. 
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GYRO 


PRECESSION $ The spin momentum of a gyro precesses 





VECTOR relative to inertial space in an attempt 

Wicge) to align itself with the applied torque.. 
A Fin 
-~ . TORQUE 
| VECTOR 

APPLIED FORCE 7 
PRODUCING TORQUE , 
, GYRO 
PRECESSION 


MOMENTUM VECTOR 
H, 


FIGURE 13-2. REPRESENTATION OF THE BASIC LAW OF MOTION 


Wigs 


Hs 


FIGURE 13-3. DEFINITION OF THE PRACTICAL GYROSCOPE AXES 
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Referring to (13-5), we express each vector by its components by 


Ww, \* | 0 
Woes Wy Hy = : 
W H, 
and 
Iw, (J, P, W,)\" (J: PW, 
_|J. PW 
PH =P,, 1 Woy oll PW, ‘ ; 
IW, IP. W, vee WY 
M. \¢ 
i M, 
M 


where Jx, Jy, and Jz are the moments of inertia of gyro-element about its x, y, and z axes. (J is 
used instead of I, which is used for identity matrix.) Substituting these into (13-5), and 


performing the cross-product Wj. se X Hs, we get the following: 


pW, +0=M,. (13-7c) 


Understanding the components are expressed in ge—frame. 


Of these, only 13-7b is of interest to us, which is repeated in a modified form with original 
subscripts for Wx and Wy. 


Py W(i-geyy)= My + Hs; W (-ge)x . (13-8) 


Equation 13-8 shows that the time rate of change of the angular momentum about the y-axis 
of gyro-element is equal to the sum of the applied torque about the y-axis and the gyro-scopic 
torque created by the interaction of the gyro’s spin angular momentum H, and the angular velocity 


of the gyro-element about its x-axis. 


Some authors refer to this gyroscope torque as the “reaction” torque caused by the 
Wa-ge)x- However, it is not “reaction” in the sense of Newton’s third law because H;W -ge)x 


acts in the y-direction, while Wq_ge)x is in the x-direction, which is perpendicular to the y- 
direction. 


13-4 
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14.0 SDOF GYRO 


In an SDOF gyro, the gyro-element (ge) consisting of the rotor and its gimbal is encased 
inside a case called the gyro-unit (gu), in such a way that the ge gimbal axis is co-axial with gu 
case's longitudinal (cylindrical) axis, which is called the output axis (OA), as shown in Figure 
14-1. Thus, this type of gyro has no gimbal other than the ge-gimbal, which holds the spin-rotor. 
So, in the SDOF gyro, the gyro spin vector (called the H vector) denoting the angular momentum 
of the spin-rotor, is allowed to precess (rotate) about the output axis of the gu (case) only. ’ 






OUTPUT AXIS 
OA 
gyro - unit ee 
én GYRO ELEMENT 
GIMBAL 
ELASTIC 
SIGNAL RESTRAINT 
GENERATOR 


FIGURE 14-1. ESSENTIAL ELEMENTS OF A SDOF GYROSCOPE 


The ge frame and gu frame are defined and depicted in Figures 14-2 and 14-3. Both 
frames share the same origin and the same ¥ axis, (i.e., Y ge = Ygu). The ge frame with spin axis 


along Zg- rotates about the common axis of Y ge = Y gu making angle & from Zgy(SRA) to Zge (SA). 


14-1 
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GYRO ELEMENT GIMBAL 






ROTOR 


ROTOR DRIVE 


FIGURE 14-2, THE GYROELEMENT 


Pa 


Yge GYRO ELEMENT 
GIMBAL AXIS 


& GYRO OUTPUT ANGLE 


Z 
gu 
SRA SPIN REFERENCE AXIS 


SA SPIN AXIS 


Xge ‘ 


Xgu 
IA INPUT AXIS 


FIGURE 14-3. GYRO ELEMENT COORDINATE AXES 


Xsgy - along IA (input axis) of gyro-unit (case) 
Ygu along OA (output/precession) axis of gu (case) 


14-2 





NSWCDD/MP-95/87 


Zg, along spin reference axis (SRA) of gu (case) 


When a= 0, Zgy = Zze or SRA =SA. 


a is the angle the spin axis (SA) makes WRS the SRA, which coincides with the Z gu) as a result of 


the rotation of the gyro-element about OA. 
Now 

Wiege = Wi-gu + Weu-ge - 
Expressing (14-1) in the ge frame: 


Wwe 


ge = Cou tga + Wo ee 
Referring to Figure 14-3 and noting that for a small angle a, sina=a and cosa = 1, 
Xge = Xgy COSA + Zgy SING = Xgy + AZgy 
Yge = Ygu 
Zge = —XgySin® + ZgyCosa = —X yy + Zey . 


From the above equations, 


Be ge gu 
Xi 1 O@ X ou 
te =} 0 1 0 13 
Lae —~ 0 1 a. Zo 
Thus, ce in (14-2) is given by 
1 O «a 
Ca = 0 1 0 
—a 0 1 


We designate the components of WF. and WE®" as follows: 


ge I-gu 
W, 
Wrese =| Wy 
WwW. 


z 





(14-1) 


(14-2) 


(14-3a) 
(14-3b) 


(14-3c) 


(14-4) 


(14-5) 


(14-6) 
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Wia 
Wre =| Woa . 
Weea (14-7) 


Now, ge can rotate with respect to gu about Ypge = Ygu = OA only, creating the angle a. 
0 
We = |Pa (14-8) 


guU'ge 


Substituting (14-5), (14-6), (14-7), and (14-8) into (14-2) : 


W, = Wia + OWsRA (14-9) 
Wy = Woa + Pa (14-10) 
W, =-aW1, + Wsra- (14-11) 


The applied torque My about the OA ( Ygu = Yge axis ), consists of : 


-Ka—Elastic torque proportional to angular displacement @ and opposing o 
-CPa—Damping ( dynamic-—friction ) torque proportional to angular velocity Pa, and 
opposing Pa 


M,,—Commanded torque applied by the torque generator located on the output axis 


Thus, we may write 
My, = —Ka — CPa + Mig (14-12) 


By designer’s choice, when Ka predominates ( relative to CPa term ), the instrument 


becomes a rate gyro. If the CPa term predominates, the instrument becomes a (rate) integrating 
gyro. This is explained below: 


In Section 12.0, we derived the following equation: 
JyPW, — H;W,x = My (14-13) 


substituting (14-9), (14-10), and (14-12) into (14-13): 
JyP(Woa + Pa) —Hs(Wig + OWspa) = —K. — CP + Mig (14-14) 


14-4 


yO 
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which gives (expressed in the ge frame): 
JyP20, = HsWia +Mig — Kas — CPx + HyoeWora — JyPWoa (14-15) 


For the operational environment in which SDOF gyros are used, (&Wspa ) and PWoa are 


negligible and, therefore, the last two terms on the RHS of (14-15) may be ignored. It follows 
from (14-15 ): 


JyP2a = HsWia — Ka — CPa + Mig. (14-16) 


Equation (14-16) shows that, in the absence of the commanded torque Mig, the gyroscopic 
torque H,Wya causes angular acceleration P2a, which overcomes the elastic torque Ka and the 
damping torque CPa. 


14.1 RATE GYROS 


For rate gyro applications, My, = 0. From (14-16): 


J H 
ZPat+ePato=W,, (14-17) 


By design K>>C>>Jy for a rate gyro. So, for a constant input value of W1a, & settles 
rather quickly to the final value of 


a= Te Wis (14-18) 


J 


Cc 
y . 
because x << <<1 in(17). 


Since the output angle @ is proportional to the angular rate Wy, of the gyro unit (case) 
relative to the inertial space, the SDOF gyro of this type is called “rate gyro.” 


14.2 RATE INTEGRATING GYRO (RIG) 


For a RIG, the elastic constant is removed, thus making K = 0, and the damping constant 
C is made large (larger than that of the rate gyro). 


Thus, from (14-16) setting K = 0, 


or 
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J, H 1 
P| Pata |=—*W, +=M,,. 
Cc C C (14-20) 


J | 
For a constant input value of Wya, the GPa term vanishes rather quickly (relative to a) 


J 
because C>>Jy by design or << 1, and (14-20) reduces to (in the absence of Mg), 


C 
H, 
poe A (14-21) 
| d | 
Since P (-) =a (-), we have from (14-21) 
Het 
i Chute (14-22) 


According to (14-22), the angular displacement @ of ge relative to gu (case) about the 
output axis is proportional to the integration with respect to time of the angular rate W,,. Hence, 
the SDOF gyro of this type is called a “rate integrating gyro” or simply an “integrating gyro.” 


t 


Denoting i} W,, dt=W,, t=0,, in (14-22) for a constant Wja, we have 
0 


a= a ~ Fig 
ee ee (14-23) 


H 
It is to be noted that Oo, is not exactly equal to O,, but scaled by Cc: 


14.2.1 Gyro Transfer Function 
(Those who are not familiar with the Laplace Transform method may skip this section) 


By Laplace Transform, differential equations are transformed to algebraic equations. This 
sometimes facilitates the system analysis, without computer simulations of the differential equation 
in the time domain. 

Denoting the Laplace Transform L of a(t) by A(s), thatis, L[a(t)] = A(s) and recalling 
that L{Pa(t)] = SA(s) assuming o (0*)=0 


and 


14-6 
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L[P2c(t)] = s2A(s) assuming 0.(0*)=0 Pa(0*)=0, we get from (14-16), with Mi = 0, 


Jys2ar(s) + Csa(s) + Ka(s) = H;Wia(s) = (JyS2 + Cs + K)a(s) = HysWya(s) (14-24) 
which gives 
( fe W,, (s) (14-25) 
a(s)=——_——- : : 
: J ; s7+Cs¢+K “ . 
or 
a(s) = G(s)Wya(s) | (14-26) 
where 
a H 
G(s)= on (S) : (14-27) 





= 


cs) J s°+Cs+K 


G(s) is called the transfer function, which relates the input Wya(s) to the output a(s) = aoa(s) 
with subscript OA indicating that o is the rotation about the output axis. 


For a constant step input of Wj,, W,, (s)= > Wa: 
Applying the final value theorem of lim  (t)= lim sa (s), we have, from (14-25), 
t-o s-0 
H, 
a (t)=lims —-————_ W,, (14-28) 


s 2 
ae ep J ys +Cs+K 


or 





a(t)= A, as t —> 00 
kK - (14-29) 
which is identical with (14-18) for the rate gyro. Recalling that K>>C>>J, for a rate gyro (that is, 


C is negligible relative to K, yet C/Jy is big), we can see from (14-28) that a(t) approaches the 
final value faster as S approaches zero than the case when C is comparable to K or bigger than K. 


For the integrating gyro K = 0, while keeping C>>J,. For this case, from (14-25) 


a.(s) Wia(s) . 


— _H, 
S(JS+C) (14-30) 


14-7 
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The transfer function in (14-30) has a pole at the origin (s = 0), and therefore, we cannot 
apply the final value theorem. 


For a constant step input of W,, ,W,, (s) == 


if H Wy 
_ Wats S aera 


W,,- So, from (14-30) 


(14-31) 
where 


Hy Wu 


oi (IS+C) S (14-32) 





By denoting the inverse Laplace Transform of F(s) by L-1F(s) = f(t) 


H, 


eel 
s-{-£)[s-o (14-33) 


It follows, using the mathematical table for Laplace Transform : 


£(t) = as ~c{eo( 1} = Hea 1-exp(-)}. 


J (14-34) 


f(t)=L" 





In a RIG C>>J by design or >> 1,s0 we may regard exp (-F to approach zero very 


quickly. It follows from (14-34) for a step input of Wy,: 
an H,Wia 

= (14-35) 

Using a Laplace Transform identity 1n (14-31) with (14-35), 


a(t) =Lo(s)=17 <F(s) 


= [ f(x)ac 


14-8 
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t 








jae H, f H, 
~ C dT=— Wy . dT= C Ww,t 
or 
H, 
a(t) = 84 (14-36) 


where 6,, is the angular rotation about IA caused by W,, for the duration of t, or 6,4 = Wi, t. 
Note that (14-36) is identical with (14-23), for the RIG. 
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13.0 STABLE PLATFORM 
(IMU PLATFORM) 


The stable platform of an LRBM is a platform, which is kept aligned with the earth- 
centered I-frames by a gimballed servo-mechanical control system, which uses gyros to detect 
rotational motion of the stable platform. Three mutually orthogonal accelerometers are rigidly 
attached to the stable platform and are used to determine the position and velocity of the missile. 
The stable platform with the accelerometers and the associated electronics is called an ae 
measurement unit IMU and is the central part of the missile guidance system. 


First, we explain the single-axis stable platform before discussing the more general three- 
axes stable platform. Consider a conceptual single-axis stable platform (plate) with the body-fixed 
axes as shown in Figure 15-1. 


Yp 


ZB 


XR, XY 


FIGURE 15-1. CONCEPTUAL SINGLE-AXIS STABLE PLATFORM 


The X-axis is constrained to be parallel to thé Xj axis of the I-frame. The plate, which 
is perpendicular to the X-axis, is subject to the rotational disturbances about the Xp-axis caused 
by unwanted disturbances. This means that, although the plate itself has a fixed orientation in the 
inertial space, the Yp-axis and Z-axis fixed on the plate may deviate from the original Yq - axis 
and Zj- axis. 


15-1 
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The control problem is how to counterbalance the deviation of the Yg-axis and Z-axis 
caused by the unwanted disturbance rotation about the Xp-axis, so that the Yp-axis and Z-axis 
return to the directions parallel to the Y-axis and Zj - axis. 


This may be achieved by a scheme as shown in Figure 15-2. The center piece of the 
scheme is an SDOF gyro with the gyro unit's IA designated as the X-axis, the OA as the Y-axis, 
and the SRA as the Z-axis. SRA is the direction of the spin axis (SA) when SA is orthogonal to 
IA. SRA is the reference direction of SA. SA is always orthogonal to OA. 





GYRO 
ELEMENT 


FIGURE 15-2. MECHANICAL SCHEME FOR SINGLE-AXIS STABLE PLATFORM 


The type of gyro used for this purpose is called a RIG or simply an integrating gyro for the 
following reasons: 


The IA of the gyro unit is constrained by the motor-shaft structure not to be rotatable with 
respect to the base except by the servo-motor control. When the base with the gyro on it (with the 


servo control off) rotates about the +X axis (IA) by an angle 0,, the SA will precess (turn) 
downward toward the +X direction, following the gyroscopic principle that the spin vector will try 
to align itself with the torque vector (the direction of rotation). 


This precession of SA rotates the gyro-element’s SA axis by an angle 0, (from SRA 
direction) about the OA (+Y direction). 


15-2 
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For the RIG, by special design construction, Q, is related to 0, with close approximation 
for small angles by: (see a previous section on SDFG) 


o, = 6, 


YC (5-1) 


in which H, is the constant spin angular momentum of the gyro spin rotor, and C is the dynamic 
friction (so called damping) coefficient of the gyro element about the output axis (OA). Equation 


15-1 shows that the output angle 0, aboutOA is proportional to @, or to the integration of the 
angular rate 0. about JA which is the input angular rate. For this reason, the device is called a 
RIG. : 


The @, about OA caused by 6, about JA is sensed by the signal generator (SG) 
mounted on the output-axis. This electrical signal from SG drives the servo-motor whose shaft 
coincides with the JA of the gyro, so that the gyro unit is rotated in the opposite direction from the 


initial 6, rotation. 


This operation turns the spin axis (SA of gyro) upward until 6, becomes zero or 6, is 
rebalanced, at which time the signal from SG becomes zero, causing the motor operation about 
IA to cease. With this control scheme, Yg and Zg axis misoriented by the initial rotation 6, 


about IA of the platform is returned to the original Y; and Zy axes by the counter-rotation (-8,) 
about the (+A) direction of the platform. 


Referring to Figure 14.2, suppose the SA drifts downwards generating @, in the absence 
of rotation about the input axis (Xp-axis), even though Yg and Zp axis point to the correct 
directions, and thus do not need correction. 


Nonetheless, in the actual system, the SG will detect 8, and make the servo-motor turn the 


platform about JA in such a way as to null the 6,. This would cause the Yg and Zp axes to 
deviate from the correct directions when correction is not necessary. This results in the platform 
drift-error caused by undesired drift of the gyro spin axis. This phenomenon is called gyro drift. 
By custom, the angular rate of gyro drift is simply called gyrodrift. 


Consider a hypothetical platform (block) with three single-axis platform control systems 
with each SDOF gyro's JA , OA, and SRA oriented as shown in Figure 15.3. 


15-3 
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Z 

Z - gyro 
IA 
SRA 
Y - gyro 
OA 
OA 
X - gyro an ¥ 
SRA 
SRA 
OA 
IA 


FIGURE 15-3. THREE-AXIS STABLE PLATFORM 


This is a hypothetical arrangement for tutorial purposes. 


Att=0, we assume the Xs, Ys, Zg axes of the platform are parallel to the X;, Y;, Z, axes of 
an inertially nonrotating frame. 


The input axis of the X-gyro designated by IA, points to the inertially fixed direction if 
there are no rotations about IA, (IA of Y-gyro) and IA, (IA of Z-gyro). Similarly, IAy points to 
the inertially fixed directions if there are no rotations about IA, and IA,, and IA, points to the 
inertially fixed directions if there are no rotations about IA, and IAy. Thus, if, at t=0 when 
IA,, IAy, and IA, are parallel to X;, Yj, and Z; axes, and if all three control systems function 
propertly, the platform would be inertially nonrotating about each of three axis and the platform 
becomes an inertially stablized platform, even under unwanted disturbances. 


The actual three-axis stable platform control system is much more complicated. For those 
who are interested in the further study are referred to the following: 


Inertial Navigation Systems 


by: Charles Broxmeyer 
McGraw-Hill, 1964 
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Gyroscopic Theory, Design and Instrumentation 
by: Wrigley, Hollister and Denhard 
The MIT Press, 1969 


Control System Design Analysis of Three-Axis Dynamic Rate Table 
by: Kee Soon Chun 

MIT Dept. of Aeronautics and Astronautics 

(SM Thesis), 1968 
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16.0 PENDULOUS INTEGRATING GYRO ACCELEROMETER (PIGA) 


Consider an SDOF RIG with a mass unbalance m located at a distance r along the SA from 
the center of the gyro (the intersection of OA and SA) as shown in Figure 16-1. 





+Y 


ae 
OUTPUT 
AXIS 


GYRO 
ELEMENT 


GYRO UNIT (GU) 


FIGURE 16-1. SDOF RIG WITH A MASS UNBALANCE m 


The motor-shaft axis coincides with the IA of the gu. The motor is fixed to the base of an 
inertially nonrotating platform. 


Suppose the platform, such as an IM of the missile, is under acceleration whose 
nongravitational portion (specific force fg) is along the -X direction. This causes the mass m to 
experience the inertial reaction force toward +X direction, generating a torque about OA, that is, 
about the +Y axis. If the gyro unit were free to rotate about the gyro IA, the spin vector H would 
rotate (precess) in the Y - Z plane about the -X axis (W},), trying to align the H-vector with the 
applied torque (vector) about the +Y-axis. The precession of H-vector is caused by the inertial 
reaction force on m, according to the gyroscopic principle. 
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However, GU (whose input axis coincides with the motor shaft) is prevented from rotating 
about its IA by the fixed motor structure. Therefore, the SA with the mass-unbalance on it has no 
choice except to move downwards, causing the SA to rotate about OA (+Y axis) creating an angle 
A from the original null position. This rotation is sensed by SG and drives the motor shaft 
through the servo motor control in such a way that the angle A is nulled (rebalanced)—in this case 


by rotating the motor shaft in the -X direction. When A reduces to zero, the output of the SG 
correspondingly reduces to zero as well, and the motor shaft rotation stops. 


For the reasons explained in the previous sections on accelerometers, fg (gravitational 
specific force) does not contribute to the generation of torque. Only fxg causes the reaction- 
specific force fp, which generates the torque. The torque generated by the inertial reaction force 
because of fyg onthe mass m located at r onthe SAis mrfp = -mrfng with the negative sign 
indicating that fp is in the opposite direction from fyg. 


From Equation (16-19) from the Section 13.0 on the SDOF gyro, 
JoaP2Ao, + CpAoa = HsWia + mrfxe(1) 

In (16-1), the applied torque mrfyg is nulled (rebalanced) by the opposing gyroscopic 
torque HsWy, caused by the rotation of GU about the negative input axis (about -X axis) by the 


servo motor, thus causing the direction of the torque generated by HsWja to point in the -Y 
direction, which is opposite the direction (+Y direction) of the torque generated by mrfyg. 





Under steady-state conditions, both PA,, == Ag A,, and Pp? Aon = A,, are zero. It 
follows from (16-1): 
H;Wia = —mrfne 

Writing W,, = = fig= al and integrating both sides from 0 to T, the computation 


period, we have 





‘That is, the angular rotation 8 about IA accumulated during time T is proportional (scaled 
by 7 ) to the velocity accumulated (during T ) as a result of the nongravitational specific force 
fyG applied in the -X direction. 
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This device senses and measures that part of the velocity, which is the integration of 
specific force (acceleration) (caused by the nongravitational specific force fyg) by means of a gyro 
with a_pendulous mass deliberately placed on the SA. Hence, it is called a pendulous integrating 
gyro accelerometer or PIGA. 


As we have seen before in the section on SDOF RIG (Equation 16-23 of Section 14.0), the 
Aoa of GE from GU is related to Oy, by | 
H 


Aon = Oy 


C 


The Apa is measurable by the signal generator attached to 0A. Thus, in practice, Vg may 
be determined by the summation of the electric pulses from the signal generator attached on 0A, 


rather than by the mechanical summation of 6,4. 


Unlike PIPA, PIGA does not suffer from the voltage and current magnitude variations. 
Also, PIGA is almost insensitive to the variations of the magnetic flux density and permeability. 
This explains why PIGA replaced PIPA in the LRBM. 
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17.0 EQUATIONS OF MOTION FOR INERTIAL NAVIGATION 


By Newton’s second law: 


2 
mP Rp =FetF yo (17-1) 


Since the external forces have to be either gravitational force (F,) or nongravitational tre 
(Fxg) in a mutually exclusive way, by denoting 


F, 


F 
fh and ae we may express (17-1) by: 


Pe 


: F,, 
PIRp == t+ =fot fins (17-2) 


Equation (17-2) is the equation of motion to implement an INS in the earth-centered I- 
frame. We have previously found that the onboard accelerometer senses nongravitational specific 
force fxg only, leaving out the gravitational specific force, fg. This necessitates the determination 
of f, (as a function of position) by the onboard computer, based on a mathematical model to 


implement Newton’s second law by onboard instruments and a computer in a self-contained way 
without external help. 


The block diagram for an INS implementation in an I-frame is shown in Figure 17-1 
below. 
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P.R,, velocity 





R,, position 









Fa Computer 
m 


PLATFORM 


Non - rotating 
w.r.t. 


Inertial space 





FIGURE 17-1. INS IMPLEMENTATION IN AN I-FRAME 


Now, to derive the equation of motion in the earth-fixed frame (that rotates with the earth 
WRT the I-frame with angular velocity Wy), we apply the theorem of Coriolis to Rz, between the 
earth-centered I-frame and the earth-fixed (and earth-centered) E-frame. 


Thus, 
P.Rg, = P2Rg, + Wi Xx Rg, ; (17-3) 
It follows 


P Rey =P I (P R,,) 


=P,(P R +W,, xR,, ) using (17-3) 


E Ep 
= P,(P.Re, ) + P:(Wee Rep) : (17-4) 


Referring to the right-hand side of (17-4), by applying the theorem of Coriolis to the first term: 
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P(PAR,,) = Ps (P.Re,)+ Wa X(P:Rep) 


= P2R,, + We x(P5Rz,) - - (17-5) 


Noting that Wy: is constant and therefore P,Wi =0, we have, from the second term on the right- 
hand side of (17-4) 


P,(W, x Rz,)=(P,W,)xR,, + We x(PR,,) 
= W, x (PRap) 
= Wy X(PeRep + We X Rgp) 
= We x(PpRz, )+ We X(We Rg) - (17-6) 
Substituting (17-2), (17-5) and (17-6) into (17-4): 
fs + fyg = PERE, + 2Wee X(PeRe,)+ Wee X(We XRz,) - (17-7) 
For the case in which the point P is stationary on earth, Re, is constant and, therefore, 


P,R,, =0 and therefore P, (P.Re )=PER 0. 


Ep = 
For this special case, (17-7) reduces to 
f, = Wr xX (Wi x Rs) — -fing . (17-8) 


As we have previously discussed in the section entitled The Prototype Linear 
Accelerometer, the accelerometer output measures -fyg. Thus, the left-hand side of (17-8) or fg - 
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Wr x (Wx X Rgp) is the vector quantity (both direction and magnitude), which a plum-bob on a 
simple pendulum at rest on earth senses. We call this vector quantity “gravity” and denote it by g. 
That is, 


g Afco- Wu x (We Rz,) - (17-9) 
Note that in (17-9) g=f, for Wr: =0 for the nonrotating earth. For the rotating earth, g 
is not equal to f, except at the north and south poles, where Wz x Rg, = 0. 
Substituting (17-9) into (17-7): 


fug + £= PzRe, +2Wa X (P.Rz, } : (17-10) 


Applying the theorem of Coriolis to (Pg Rep) between the E-frame and N-frame where 
the N-frame is the local-level northeast down frame, 


PRe = P,(P.R,, 
= P,(PRe,)+ Wen x (PaRep) (17-11) 


Substituting (17-11) into (17-10), and coordinating (expressing vector-components) in the N- 
frame: 


[Px (PeRe; \P = gX +68, —-(2WN+W,)x(PR.) - 17-12) 


Equation (17-12) is the equation of motion for the near-earth (land, air, and sea) navigation 
implementing a local-level (north, east, and down (NED) platform. It implies that the velocity 
relative to the earth, P;Rg, and the position relative to the earth, R,, (which is the integration of 


PeRgp ) may be determined based on the f. which is the negatives of the outputs of three 
accelerometers along the NED directions on the platform, which implement the N-frame. 


Referring to left-hand side (LHS) of (17-12), we see that 
[Pu(PERe \P=Py (P.Re \" 
=P(P.P.). 
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(), 


e/a. 


Therefore, since p(-)= 


JP Re "a= JF Rea 


=(P,Ryp \"x Vv, (17-13) 


where Vx, Vy, and V, are respectively NED velocities of the vehicle relative to the earth-fixed 
frame. 


The vehicle's position relative to the earth (latitude and longitude) is determined from the 
velocity output by integration in the following way. 


For the spherical earth, denoting latitude by and longitude by A, and the magnitude of 
Rep by R, we get (refer to Figure 17-2): 


Vn 
Vi=OR or d= so that 


f ar=f oa 





(17-14) 
=o 
= latitude 
dV.=(R A or X “s 
and V,,=(Rceoso )A or ra 
so that 
V, | 
[seat = | jdt = 2. = longitude. (17-15) 
Rcosd 
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(East into paper) 


FIGURE 17-2. DETERMINATION OF LINEAR VELOCITIES FROM THE 
LATITUDE AND LONGITUDE RATES 


Next, by denoting the magnitude of Wy by Q, we have (refer to Figure 17-3): 


We Qcosd 
We =| We, [=| 0 |. 
We: —Qsino (17-16) 


The negative sign in Wy is necessary because the vertical component of W;, is upward 
and, therefore, opposite to the Z - direction, which is downward. 
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FIGURE 17-3. DECOMPOSITION OF THE EARTH RATE INTO NED COMPONENTS 
, For Wi (angular velocity of the local-level NED frame relative to the earth-fixed frame), 


the N-frame rotates relative to the E-frame caused either by latitude rate @ or longitude rate 1. 
) Referring to Figure 17-4, the North component is Acoso and the Z (down) component is — Asind , 

because Asin is in the upward direction, which is opposite the positive Z (downward) direction. 
The north and down components are not influenced by @ because @ is perpendicular to the 


| meridian plane. The east direction is opposite to the @ direction (west). Therefore, the east 
component — @ is not influenced by A because @ is perpendicular to’. Thus, 


Wenx Acosd 
) Was =| Wa |=] ~4 
, Wenz) \—Asing (17-17) 
| 
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FIGURE 17-4. DECOMPOSITION OF THE LONGITUDE RATE INTO NORTH AND DOWN COMPONENTS 


To maintain the platform (on which NED accelerometers are mounted) in the local-level 
(NED) frame, it has to be torqued by the angular rate Wi, ; 


Using (17-16) and (17-17): 


Wa = Wa t+ Wax 


Qcosd Acoso 
0 |+) -o 
—Qsin o ~r sin > 
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fe + i) cos 


= — : ‘A 
(0+4)sing (17-18) 


The implementation of Equation (17-12) for terrestrial navigation using the local-level 
(NED) frame is shown in Figures 17-5 and 17-6. 











> Latitude 
A. Longitude 
d Depth 


Gravity 
Computer 


; Latitude ¢ > 
Vv Vv A 
U Longinde A 
; | Computer : 
Platform . 
Coriolis 
Pee fC“ 
Computer 
Earth Rate 
Physical Signals Gyro - Torque 
Computer Barth Rate 


FIGURE 17-5. SIGNAL FLOW BLOCK DIAGRAM OF LOCAL-LEVEL N-FRAME 
IMPLEMENTATION OF TERRESTRIAL NAVIGATION 
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Latitude d 
Longitude A 
Depth’ d 
(Height -d) 
Coriolis Force Computer 
(2WE + Wax) x (PaRee iy > 





[W,,|= Q = constant 


Wx 
Computer 


FIGURE 17-6. IMPLEMENTATION OF NAVIGATION EQUATION (17-12) 


Repeating Equation (17-12) : 


[Ps (PaRs, )} = g" + fro — (2Wa + Wax) s (P.R,,) (17-12) 
where 
Ne Vic 
(PLRe \"= Vv, =| Vz 
v. Vp 
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| V, Vv. 
| f, g 
f" =| f, g =! 8, 
} f, g. 
| Wi Wo, 
Wa i Way Wan a7 Waxy 
Was Wa 
then 
(2W% + WX.) x(P,R,,) 
i j k 
=|2Wae, + Win, 2Wiey + Wen, 2We, + Wen] i, j, k of N - frame 
V, V, V, 
(2W,, + Wan, )V. —(2Wa, + Wen) V, 
=| (2W,, + Wa, )V, —(2W., + Wen, )V, 
(2Wre + Ware) V, —(2Wy, + Wen, )V, 
substituting these into (17-12): 
V, =f, + 8, +(2Wy, + Wan.) V, —(2Wa, + Wan, )V, 


V, =f, +g, +(2We, + Wax, )V, -(2W,, + Win, )V, 
17-11 
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ta 


in which, referring to (17-16) and (17-17), 


Wa. = S2cosh 


W,, =0 


WR = ~Qsingd 


Wax = Acos 


17-12 


V, =f, +g, +(2Wa, + Waxy )V, —(2Wax + Wen) Vy 
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18.0 DETERMINATION OF INITIAL POSITION AND INITIAL VELOCITY 


In inertial navigation, we determine the acceleration P? R,, from 


P?Ryp(t) = fg + tne (18-1) 
where 


Rep is the displacement vector from the center of the earth to the point of accelerometer 
location, and fg is obtained from the mathematical model as a function of position and fxg from the 


accelerometer measurements. The velocity P;Rep is obtained by integrating Pe Rp that is, 


PRyp(t) = [ PPRep(t)dt + PRyp (0) (18-2) 


where P;Rep(0) is P;Rgp at t= 0. The position Rgp is obtained by integrating P,Rzp that is, 


Rwo(t) = J, PRe»(B)dB + Rup(0) (18-3) 
where Regp(0) is Rep at t = 0. 
Thus, we see that, to determine the current position Rgp(t), we need the initial velocity 
P;Rep(0) and the initial position Rgp(0). 
18.1 DETERMINATION OF THE INITIAL POSITION 


We start with 


Ryp = Rus + Rgp (18-4) 
where 


E = Center of the earth 
S = Center of the submarine (ship or aircraft) INS 
P = Center of the missile guidance system. 


and the superscript I indicates that all components are expressed in the I-frame. In (18-4), the 


submarine's location Rgs is either available as an output from its navigation system or may be 
computed from it. 
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The position Re, of the missile guidance system, before launch, while in the submarine or 


aircraft is given in the ship-fixed S-frame. The origin of the S-frame is at the center of the ship 
navigation system. The x-axis of the frame is along the longitudinal axis on the plane paralleled to 
the main (or conceptual) deck through the origin of the frame. The y-axis is perpendicular, on the 
plane, to the x-axis on the ere The z-axis is perpendicular to both x and y axes ACCOrcine to the 


right-hand rule. To express Rep as Ry, , we have to know Cc. so that, 


R..=C.Re (18-5) 


The coordinate transformation matrix C from the S-frame to the I-frame is determined 
from the gimbal angles of the submarine (or aircraft) INS (SINS) if it implements the I-frame. If it 
implements the navigation frame (NED), Go is first determined from the gimbal angles of SINS, 


and then C,, is obtained by 


Cy = C, CE. (18-6) 


The coordinate transformation matrix C’, from the N-frame to I-frame is given in terms of 
N £1 


the submarine's latitude, @ (an output from the SINS), and the celestrial longitude, Ac, by (with the 
derivation to follow): 


—sindcosA, —sinA, —cosdcosA, 


Cl=|-—singsinA, cosA, —cosdsindA, (18-7) 
cos 0 —sind 
where 
Ac=Ag +A 
in which 


Ag = Celestial longitude of Greenwich meridian 
X. = Terrestrial longitude of the submarine (relative to the Greenwich meridian) 
Ac = Celestial longitude of the submarine (relative to the inertially fixed reference) 


18.2 DERIVATION OF C,, 
This derivation is shown here in detail as a sample case of determination of coordinate 


transformation matrix for the tutorial purpose. 


It turns out that it is easier to get Cc first, and then get C,, by transposing C. or 
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To get Cy (refer to Figures 18-1a and 18-1b), rotate I-frame (with Xj, Yi, Zy axes) about Z; axis 


by an angle of A, (celestial longitude). Let us call the resulting frame the A-frame (with X;, Y3, 


Zn axes). 


+ 


This operation yields (from Figures 18-1a and 18-1b) : 


or 


X,Y (cost sind OY/(X,) 


Y, | =|-sinA cosA O|| Y, | 
Z, 0 0 1)\Z, (18-8) 
R*=CiR' (18-9) 

Zy =, 





Equatorial 
plane 





FIGURE 18-1A. DETERMINATION OF C* 
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X, = XycosA, + YysinA, 
X, = —XysinA, + YycosA, 
Z, = Z] 





polar axis 
out of paper) Yn 
Zi= Z), YJ 
Equatorial 
plane 
X] Xi 


FIGURE 18-1B. DETERMINATION OFC; 


Next, rotate A - frame about the negative Y, axis by an angle of 6 (latitude) as shown in 
Figure 18-2. 


Let us call the resulting frame the ® - frame (with X9, Yo and Zy axes). 


This operation yields: 
® é ® A 
X, cosod 60 sing) (xX, 
Y,| =| 0 1 0 Y, 
Z, —sind Q cosd) \ Z, (18-10) 
or 


R® =C°R* (18-11) 


18-4 


mer Ee CED, Seppe, ES, Pen a, AEE mae mera, moe 


NSWCDD/MP-95/87 


parallel to 
west axis 


(out of paper) 





FIGURE 18-2A. DETERMINATION OF C* 
Xy = X,cosM + Z,sing 
-Yo=—-Y, or Yo= Y) 


Zo =—-X,sing + Z,cosh 





-Yy, 
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Referring to Figure 18-3 below, we see that the N-frame (with Xjy (north), Yn (east), Zn (down) 


axes) is related to ® - frame (shifted from the center of the earth to a parallel frame located at the 
vehicle position on the meridian ) by 


Xn=Zy; YN=Yo; ZN=Xp 


Or 
Xv) (0 0 1\¥fx,) 
Yy} =| 90 1 OF] Y, 
Zy} \-1 0 0),(Z, (18-12) 
or 
R™=C)R® (18-13) 
XN Zp North 
North Xo 
up 
YN East 
into paper Yo East (West 
into paper out of paper) 
ZN 
down 
meridian meridian 
plane plane 


FIGURE 18-3. DETERMINATION OF Cc, 
Thus, from (18-13), (18-11), and (18-9) 
R® =C*R® 


= CXC°R" 
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=CC CR. 3 | (18-14) 
Since 
R* =C*R! 
we have 
NAN wv 
Cr c,cce. (18-15) 


It follows using (18-8), (18-10), and (18-12): 


0 0 1)fcoso 0 sing\’(cosA sind OY 
cr=|0 1 0 0 1 O —sind cosa 0 
-1 0 0O),\-sind 0 coso) \ 0 0 1) 
or 
—sindcosA, —singsind, cos)” 
CX=| -sind, cosA, 0 


—cosocosA —cosdosind, —sing), (18-16) 


Taking the transpose of (16) (exchanging rows and columns) will give Cc. as given in (18-7). 


18.3 DETERMINATION OF INITIAL VELOCITY PiREP(0) 
Applying Coriolis law to Rgp between the I-frame and the E-frame (see Section 7.0): 
P;Rep = PgRgp + Wye X Rep. (18-17) 
Substituting (18-4) into the right-hand side of (18-17) for Rep, 


P;Rep = Pg(Res + Rsp) + Wiz X (Res + Rop) 


=PERges + PeRsp + Wr X Res + Wy X Rop. (18-18) 


Referring to the second term on the right-hand side of (18-18), applying the Coriolis law between 
the E-frame and the S-frame, 
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PeRsp = PsRsp + Wes X Rsp 
= Wes X Rsp (18-19) 
because Rsp is fixed in the S-frame, thus making PsRsp = 0. 


Substituting (19 18-) in (18-18) and noting that Wr + Wes = Wjs, we have 
PyRep = PeRes + Wes X Rsp + Wie X Res + Wie X Rsp 


= PeRgs + We X Res + Ws X Rsp. (18-20) 


Evaluating (18-20) at launch (t = 0) and expressing it in the I-frame because that is the frame we 


use for the missile guidance and navigation, we have 


(PiRep )p=(PeRes p+ (Wu: * Res), + (Wis * Rep), (18-21) 


Equation (18-16) shows that the velocity of the missile at launch relative to the center of the earth in 
the inertial space may be determined in terms of the values of the three terms on the right-hand side 
of (18-21). 


We deal with them one by one. 


First, (P,R,, ), is the initial velocity of submarine or aircraft relative to the earth-fixed frame with 


components expressed in the I-frame. The term (P,R ‘’ is the initial value (at t = 0) of the 
EES /9 


submarine's velocity relative to the earth-fixed frame expressed in the N-frame. It is an output 
from INS and is transformed from the N-frame to the I-frame by 


(PaRus), = C3 (0)(P2Ras), 
where C,,(0) is the value of C,, att=0. 


Next, for (W,, xR, ), , we determine Ww. x Ri in E-frame and transform it to the I-frame. In the 
E-frame 18- 


0 
Wy, =| © | in which |W,, |= W 
‘W 
where W is the magnitude of the earth-rate. 
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Thus, 
i j —k 
(W, xRi;),=|0 0 w 
X, Y, Z, E-frame (18-24) 
= —1Wyo + JwXxo 
in which 
X, 
R®, =| Y, 
Z, 
is available at an intermediate point of the computation from INS. 
Then (W, xR,, ), is obtained simply by 
I E 
(Wa x Rys), = Cz (0)(War x Res), (18-25) 


where Cc. (0) is the coordinate transformation matrix from the E-frame to I-frame at launch (at t = 
0). 


In our case both the I-frame and the E-frame are polar-geocentric, thus sharing the Z (polar) axis, 
Thus Ci is 


E 


cOoSwt —sinwt 0 

C,(t)=| sinwt coswt 0 
0 0 1 (18-26) 
derived by rotating the E-frame about the polar axis relative to the I-frame by an angle of Wr, 
where w is the magnitude of the earth rate and t is the time elapsed since the two frames coincided 


(i.e., Xp axis with X] axis and Yg axis with Yg axis while the Zg axis always remains identical 
with the Zi axiS). 


Finally to evaluate the last term on the R.H.S of (18-21), (W,, x R,,, is , we determine it first in 
S—frame and then transform it to the I-frame via the N-frame. That is, 


(W, x R,, ) = c (O)CE (0)( Ws x R,, ) (18-27) 
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In (18-27), Rg, is simply the missile (its guidance system's) location in the ship expressed in 
S—frame. Wr is determined by the guidance system by differentiating with respect to time its 
gimbal angles, which are the measure of the angular deviations between the guidance system's 
inertially nonrotating frame and the missile structure at launch time, which is fixed to the ship. C5 (0) 


is the coordinate transformation matrix from the S-frame to the N-frame at launch, which is 
determined based on the SINs gimbal angles (which is the measure of angular deviation between 


S-frame and N-frame). C., (0) is determined from the latitude and celestial longitude as given in 
(18-7). 
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) 19.0 CROSS PRODUCT STEERING 
(FOR ROCKET VEHICLE GUIDANCE) 


Let us denote the missile’s current velocity by Vy, the required velocity (such as correlated 
velocity) by Vp, and the velocity to be gained (or to go) by Vg, so that Vjy becomes Vp as Vg 
goes to zero. | 
Then, 


Vg=Vr- Vu. | (19-1) 


Thus, if we drive Vg to zero with thrust, then Vyy approaches Vr. It is intuitively obvious 
that Vg will decrease most effectively if thrust is directed parallel to the Vg direction. 


In fact, it turns out that, in all current high-acceleration rocket missiles, if we steer (rotate) 
the missile’s thrust vector with the same angular velocity w as the Vg vector, then Vg will be 


driven to zero, thus guaranteeing that the missile will achieve the required velocity Vr. To obtain 
the mathematical expression for this w, denote the unit vector along Vg by u, 





or 
Vv, 
u= 
| V, (19-2) 
By differenting (19-3) WRT time: 
i VV -V.V, == V, Vo _ V.V, 
= Vv. V_ V. (19-4) 


Since 1 is a unit vector, its magnitude remains to be one, while its direction may change. 
Thus, by denoting the angular velocity of the u-vector with respect to the inertial space by w, we 
have, as derived in the section on theorem of Coriolis, 


u=Wwxu (19-5) 


where w is the angular velocity vector of u with its rotation axis perpendicular to u . 


It follows from (19-5): 


uxu=ux(wxu) (19-6) 


19-1 
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Using the vector identity 
V, x(V. XV.) = (V.-¥,)¥. -(¥,-¥.)W3 (19-7) 
we have from (19-6): 
uxu=(u-ujw-(u-w)u=w (19-8) 


because u-u=1 and u- w=0 since w Is perpendicular to u as mentioned above. 


Substituting (19-3) and (19-4) into (19-8): 





G Vi; Vo (19-9) 
= Vere Vea Nae Mae, 
Y.vev.. OV: (19-10) 
The second term of (19-10) is zero because V, X V, =0. It follows: 
V. x Vo 
w= A 
V? (19-11) 


In (19-5), (19-8), and (19-11), we conclude that w is proportional to the angular velocity of the 
unit vector Vc. Equation (19-11) shows that w may be determined from Vz, and V,. And, if 


the rocket’s thrust vector is steered with the angular velocity of w as determined by (19-11), Vg 
will be most effectively nulled assuring Vy will become Vr. 


In practice, (19-11) is used in the form of 
w=K,(V,xWV,) (19-12) 
where Kj is adjusted for good performance and may or may not be a function of Va. 


Denoting the missile’s roll axis by igx, pitch axis by jpy and yaw axis by kgz, in the 
missile’s body-fixed frame, we have 


19-2 
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Vo x Vo =|Vox Vey Vaz 
Vox Vey Vez 


= iy (Voy Voz — VozVor) 

+}p (Voz Vox 7 Vox Vez) (19-13) 

+k3(Vox Voy — Voy Vex) - 

Thus, the commanded angular velocity in pitch and yaw directions are: 
Wich = Ky (Vox Voz — VozVox) along jg - axis (19-14) 
Wyaw = Ki (Vex Vey — VavVox) along ky - axis (19-15) 

ie does not effect the direction of thrust, and, therefore, is not used. 

In the LRBM, w is determined by 


W =K,(A, x Va) (19-16) 


instead of (19-12), where A; is the nongravitational specific force fyg, which is equal to the thrust 
in a vacuum. To see that (19-16) is a modification of (19-12), we start with the differentiation of 
(19-1): 


V, = Vato (19-17) 
Since 
V, =gt+f, =f+A, (19-18) 


and noting that V, is a very slowly varying function and thus V, may be approximated by zero, 
and also noting that |A,|>>|g| near the deployment release point we have from (19-17) and 
(19-18): 


Vo=—-Vy =—-g-A, =—-A;.- (19-19) 
Substituting (19-19) into (19-12): 
W =K(V, x (-A;)) = K(t)(A; x Vq) (19-20) 


which is the same as (19-12), with K now being a function of time for optimum performance. 


19-3 
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20.0 DERIVATION OF WEIGHTING MATRIX W 
FOR STELLAR SIGHTING 


20.1 INTRODUCTION 


The position and velocity of the missile as determined by the onboard guidance system 
during powered flight have inherent errors caused by the errors, which exists in the accelerometers 
and gyros, as well as those caused by the launch-time initial conditions associated with navigation 
and fire control. 


By means of a stellar sighting of a predesignated star before to the initiation of the 
deployment of the reentry bodies, the guidance system measures the change in elevation AE and 
change in bearing AB of the star relative to the guidance system platform axes. 


The guidance system uses these measurements to make the best estimate, in the least square 
sense, of the errors in position, velocity, and platform orientation by using a weighting matrix 
denoted as the W matrix. The W matrix is predetermined before launch, based on the error 
Statistics of the accelerometers and gyros as well as those of the initial conditions related to 
navigation and the fire control. 


20.2 FORMULATION 


Define a state vector x, which is to be estimated at the time t, of the stellar sighting by: 


X =| Ap, (20-1) 


where Ap,, Ap, and Ap3 are the guidance system's position errors, Ap, , Ap,, and Ap, are 


velocity errors and 8}, 82, and 83 are platform tilt errors, along the X, Y, and Z axes, denoted by 
subscripts 1, 2, and 3 for convenience. 


We denote a constant error vector with some 65 elements by e. Each element of the e 
vector is considered as a random constant for a given missile flight, but varies from flight-to-flight 
randomly with known standard deviations and mean values. 
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We relate the state X at the time of sighting to e by means of the M-matrix, that is X=Me. 
Note that X being a (9 by 1) column vector and e being (65 by 1) column vector, the size of the M- 
matrix must be (9 by 65) to be conformable. 


In the stellar sighting, the onboard instrument measures the elevation error AE, and the 


bearing error AB (of the designated star) relative to the guidance system's platform axes. We 
represent this by the 2x1 measurement vector Z by: 


( 
—- 
AB (20-2) 


Z=HX+v (20-3) 


We relate Z to the state vector X by: 


where v is a 2x1 measurement noise vector, and H is called the measurement matrix. Since Zis a 
2x1 matrix and X is a 9x1 matrix, H must be a 2x9 matrix to be conformable. 


Our objective is to find the weighting matrix W such that the estimate X of X based 


on the measurement Z is optimal (in a sense to be defined later). That is, we have to determine the 
weighting matrix W such that X in 


X= WZ (20-4) 
is optimal. 


It follows, substituting (20-4) into (20-5): 


X = W(HX +v) (20-5) 


20.3 DETERMINATION OF W MATRIX 


The state vector X represents the true value of the errors in the guidance system's 
determination of position, velocity, and platform alignment. We seek the best (optimal) estimate X 


of these guidance system errors. Denoting by X the error in the process of obtaining the estimate X 
of the system error, we may write 


X=X+X (20-6) 


It follows using (2) and (5) in (6) noting X = Me: 


20-2 
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K=XK-X 

= W(HX + v) - X 

=WHMe + Wv - Me 

= (WH - I) Me + Wv (20-7) 
where I is an identitiy matrix. We denote the covariance matrix of X by Px . That is: 

P; AE(XX") (20-8) 

Where E denotes the statistical expectation and the superscript T denotes a transpose. 
Substituting (20-7) into (20-8): 


Py = E(XX7) 


{[(WH — 1)Me + Wv][(WH - I)Me + Wv]'} 


E 
E{[(WH - I)Me + Wv][e™M"(WH-1)? +v"w" |! (20-9) 
E 


{[(WH - I)Mee"M7(WH-1)"]+ Www" } 


with the assumption that 


E(ve") 0 
E(ev") = ; 


This is a reasonable assumption because in reality, e and v are unlikely to be correlated. 
Note that the expectation E operates only on the random vectors such as v and e, but not on the 
constant matrices such as W, M, or (WH-I). 


(20-10) 


From (20-10): 
Py = E[(WH —I)ME(ee™)M"(WH - I)" + WE(vv")W"| 


= (WH -I)MP,M"(WH-I)' + WP, W" (20-11) 


where 
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P, AE(ee") 


(20-12) 
Py AE(vv" ) 


We want to determine W such that the scalar sum of the diagonal elements of P, called trace (tr) of 
square matrix P, is minimized because, generally, the diagonal elements are most significant and 
predominant. To achieve this, we set | 


iP _ 0 (20-13) 


and we solve the resulting equation for W. 


For those who are interested in a tutorial derivation of (20-16) to (20-22), see the last paragraph of 
this report. 


From (20-11) 
Py = (WH~1I)(MP,M*H™W" — MP,M")+ WP, W, 


= W(HMP,M"H")W" —-MP,M"H™W" 
—~WHMP,M‘ + MP,M’ -- WP, W* (20-14) 
Now, since P, is a covariance matrix, P, =P) and we have: | 
(HMP,M™H")" =(M"H™)’ (HMP,)" 
= (H") (M") P:(HM") 
= HMP,M"H" (20-15) 
which shows it is a symmetric matrix. 


According to a matrix formula , denoting trace by tr: 


0 


say t (WBW")=W(B+B") 


=2WB if B=B' (20-16) 


Thus: 


20-4 
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of W(HMP,M"H" )W? = 2W(HMP;M‘H") | 


ow (20-17) 
and 
d z (20-18) 
swt WP, W =2WP, | 
‘ T 
since P,=P,. 
Now, by formula: 
ae AW'=A (20- is 
ow 
Thus 
0 TryT T TryT 
att (MP_M"H")W" = MP.M"H (20-20) 
and also, by formula: 
~~ ewe=ct (20-21) 
OW 
Thus: 
—-¢ W (HMP,M") = (HMP,M")" 
ow E E 
=(M")'(HMP,)" 
=M(P,)' (HM)" 
= MP,M"H™ (20-22) 
since P, = Pi. 
Note that, since MP, M’ is not a function of W, 
d T 
=—tr(MP,M")=0. (20-23) 


oW 
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Thus, from (20-14) using the above equations, and setting the result equal to zero: 


a oP, = 2W,HMP,M'H? —2MP,M'H" +2W,P, =0 (20-24) 


with W replaced by Wo, the optimal value of W. 


Solving (20-24) for W,: 


W,(HMP;M"H" + P,) = MP;M"H" (20-25) 

Or 
W, = MP.M™H"(HMP;M‘H" +P,) (20-26) 

or 
W, = MP; (HM)"| (HM), (HM)' +P, "] (20-27) 


if the inverse exists. This is the W-matrix we are seeking. Further analysis shows that Wo given 
in (20-26) is indeed a minimum as required. 


The optimal estimate X is obtained by (20-5): 
X= W,Z 
with W, given by (20-26) 
For the tutorial, heuristic derivations of the matrix formula used in this section, 
demonstrated by means of 2x2 matrices for satisfaction (not as a rigorous derivation), the readers 
may refer to document NAVSWC MP 91-03 "Derivation of Discrete Kalman Filter Equations 


using Heuristic, Tutorial Approach" (Revision A) by Kee Soon Chun, NSWC, Jan 1991. Anyone 
interested in having a copy may contact the author at NSWCDD (Code K44). 


20-6 
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21.0 SCHULER PERIOD AND SCHULER TUNING 


Consider a spherical, nonrotating earth. If one stands still on the earth and holds a plumb- 
bob, it will point to the center of the earth in the vertical direction. If one accelerates, we know by 
experience, the line of the plumb-bob will lag behind (from the stationary direction) because of the 
inertial reaction. Thus, we can see that the plane perpendicular to the plumb-bob line cannot serve 
as the true horizontal plane because when accelerating, the plumb-bob line deviates from the true 
vertical. That is, the reference frame constructed based on the plumb-bob is not accurate because 
the frame tilts under acceleration. i 


Question : How can we make the plumb-bob follow the vertical even under acceleration? 


It is intuitively obvious that, if we can make the angular rotation of the plumb-bob line 
relative to its pivot point equal to the angular rotation of the earth’s radius vector from the center of 
the earth to the pivot point, the problem is solved. 


Now, consider a hypothetical simple pendulum with the length equal to the radius of the 
earth and the mass m located at the center of the earth. It is obvious that, no matter how the pivot 
accelerates, the line of the pendulum always points to the center of the earth along the line of 
vertical. 


From elementary physics, we know the period T of a simple pendulum with length L is: 


T=2n/= 
£ 


For the pendulum with L=R, the radius of the earth: 


T=2n | = 84min 
& 


In his classical paper (published in 1923) on the research of the gyrocompass (North- 
seeking device), a German scientist, Max Schuler suggested that if we could build a servo- 
mechanical control device (platform), which has a natural period of oscillation of 84 min, the 
device will track the earth’s vertical under acceleration without oscillation, if it was initially aligned 
with the earth’s vertical. 


For this reason, the 84-min period is referred to as the Schuler period, and the natural 
frequency corresponding to the Schuler period is called the Schuler frequency. That is, 


W, =2nt=2n— = 
T 2m VR 
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Ww. = 1 = is the Schular frequency. 


A physical device (such as an inertial platform) that has this arrangement is said to be 
“Schuler tuned.” | 

If the platform is Schuler tuned this way, the angular velocity of the platform’s vertical axis 
is equal to the angular velocity of the true vertical, if it is intially aligned. If it is not initially 
aligned, the platform’s vertical axis will oscillate about the true vertical with a natural frequency of 
an 84-min. period and with an amplitude of oscillation equal to the initial angle of misalignment. 


21-2 


——————— 
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22.0 CORRELATED VELOCITY 


22.1 INTRODUCTION 


The correlated velocity is the required velocity of a reentry vehicle (RV), which will hit the 
target on the rotating earth by free-falling in a specified time of flight (tg). It is the required velocity 
correlated with the target displacement in the inertial space caused by the rotation of earth during 
the specified TOF of the RV. : 


The specification of tg is necessary because the target is located on the rotating earth and 
thus moves to a new location in the inertial space during ts, and because we want the RV to hit the 
same target at the new location in the inertial space. 


Another way of looking at the problem is how much speed and in what direction the RV 
should have so that, by free-falling along a elliptic trajectory, it will intercept the target, which is 
moving along a circular orbit in the inertial space, while still fixed on the earth. 


Referring to Figure 22-1, the problem of finding the correlated velocity may be stated as 
follows: 


Given: r(0) =1o at t=0 of RV, 
ry(O) of target at t=0 fixed on the rotating earth 
te time of flight of free-falling RV released at t=O until it hit the target 


Find: V, (Correlated Velocity magnitude) and y, (flight path angle) such that 
I (te) =Lr at t= te 
that is, such that the position of RV coincide with the position of target on the 
rotating earth at t= te or, RV intercepts the target 


This problem is known as Lambert’s problem or the Lambert Guidance problem. 
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FIGURE 22-1. CORRELATED VELOCITY 


The derivation of equation to determine the correlated velocity, in the approach adopted in 
this report, in terms of ro and rp(ts) is based on two key equations. 


The first equation is the expression for the conservation of angular momentum along the 
line normal to the plane containing the RV’s trajectory. The second equation is the equation of 
motion for the RV derivable from Newton’s second law of motion along the radial direction of the 


RV from the center of the earth. 


22.2 CONSERVATION OF ANGULAR MOMENTUM (refer to Section 1.0 for notations) 


Applying Newton’s second law to the RV mass, m 





r 
es (22-1) 


22-2 


aes eee eee 


sdtmcnemectntcel (Bienen fmocnimmais * TG, erin aise ee Gas, acai, 
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where P,(-) represents < (-) in inertial space and M is the mass of the earth that is assumed to 


be concentrated at its center. Using P r=i, Pr r=f and ut = GM in (22-1): 


Pr=¢=- (22-2) 
r 


multiplying both sides of (22-2) by r x (vector cross-product): 


rxt=-Gr (22-3) 
Noting that rxr =0, we have from (22-3): 
rxrf=0 (22-4) 
Since, by simple differentiation 
de 3 
qrixtaixitrxe (22-5) 
and since tx r=Q, we have using (22-4): 
d 
aq (Ext)=0 (22-6) 
or 
rt xr=h=constant vector. (22-7) 
This means for anytime t, 
£(0) x£(0)=r(t) xt (t) (22-8) 


To proceed further, we define a new rotating frame of reference that we call the T-frame. 
Referring to Figure 22-1, the x-axis is along the r(t) direction and the y-axis is along the line 
perpendicular (in counter-clockwise direction) to r(t) on the trajectory plane. The z-axis is normal 
to the x-y plane (trajectory plane) with the origin at the center of the earth sharing the origin with 
the earth-centered I-frame. The T-frame rotates about the z-axis with the angular velocity W,.= 6 


of the T-frame relative to the earth centered I-frame. 


‘Since the h vector is by definition of cross product, perpendicular to both randf, itis 
perpendicular to the plane formed by randi. Since h is constant, the plane containing randi 
is constant as well for all t, which means that the trajectory of RV remain on a fixed plane. 
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Applying the law of Coriolis to r(t) between the I-frame and the T-frame (refer to the 
section on the law of Coriolis for derivation): 


Pir(t)=P r(t)+W,xr(t) — (22-9) 


where 


d 
i = (-)— time rate of change of vector (-) as observed in the I-frame. 


.. 


P.(‘)= a ( -)— time rate of change of vector (-) as observed in the T-frame. 


In the T-frame, according to our definition: (Note: superscript T should not be confused with 
transpose.) 


r \? t \t 
r(t)=|9] ; Pir(t)=|9 (22-10) 
0 0 
0 T 
T 
Wis . (22-11) 
so that 
ijk 0\F 
W,,x2r(t)=/0 0 6) = |r6| expressed in T-frame . (22-12) 
rdQq T 0 


22-4 
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Substituting (22-10) and (22-12) into (22-9) : 
t T 0 T r T 
[Pr]'=2'=|9| +|r6 | =|r6 (22-13) 
97 \0) \0 


with components expressed in the trajectory frame (T-frame). 


It follows, expressed in the T-frame: 


ij k QO \F 
h’=r™zt=it O O} —| 0 (22-14) 
a r r6 Olr r’6 


According to (22-14), the angular momentum viewed from the trajectory plane is pointed 
along the T-frame’s z-axis as it should be and is equal to 


Daal 2 5 — 2 dO 
h=r°0=r qe (22-15) 
From (22-15), we get 
rr 
dt=+- dé (22-16) 


(22-16) is the first of two key equations necessary to derive the equation to determine the correlated 
velocity in terms of known quantities. 


EQUATION OF MOTION—along the x-axis (radial direction) of the trajectory (T) frame. 
Operating on both sides of (22-9) by P,: 


P, (Pr\=P, (P,r)+P,(W, xr). (22-17) 


Referring to the right-hand side of (22-17), treating Pr and Wy, xr as vectors and applying 
Coriolis law: 
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P (P,r)= Pp. (P.r) + Wr x (Pr) 
= Port+W,, x (Pr) 


and 


Pp 


I 


(Wixr)=P, (W,,xr)+W, x (W,, xr) 
=(P,W,,)xr+W x (P r 


Ir T 


Substituting (22-18) and (22-19) into (22-17): 


Pir=Pirt (P,W,,)xr+2W,,x (P,r)+ Wx (Wr). 


0 \t ty 
Noting that 6 is along the z-axis and thus W,,= 0 , and r= o , we have 
0 


ijk) (0, 
(P,W,,)xr=|0 0 6) ~|r6 
r 0 Olr 0 


ijk| (oy; 
W,,*(P,r)=/0 0 6 -(s 
r00l, (0 


and using (22-12) 


22-6 


)+W,,x (W,, xT). 





(22-18) 


(22-19) 


(22-20) 


(22-21) 


(22-22) 


(22-23) 


(22-24) 





NSWCDD/MP-95/87 


i j k _r62\" 
W,,* (Wy xT )= 00 6] =] o 
0716 0|, \ 0 


Substituting (22-21), (22-23), (22-24), and (22-25) into (22-20): 


t\ (0 0 —r6” 
(P?r)'= O j+J7r6 1+ 1276 [+] 0 
0} \Q 0 0 


(22-25) 


(22-26) 


with components expressed in the T-frame. It follows that the x-component (in the radial direction 


of (22-26) is: 
2_\T 42 
(P’ r } =f—r6 
and the y-component is 


(P?r \'=16 + 270 


Since the gravitational force is along the radial or x-axis: 


r. : Ss 
where 2 is a unit vector along r direction. 


From (22-2), (22-27), and (22-29): 


22-7 





(22-27) 


(22-28) 


(22-29) 


(22-30) 
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in the radial direction or along the x-axis. 
And, from (22-2), (22-28), and (22-29): 

16 + 276=0 (22-31) 
in the direction perpendicular to the radial direction, or along the y-axis direction. 


Equation (22-30) is the equation of motion of the RV along the radial direction, the second 
key equation we are looking for. 


Referring to (22-31): 
r6 +276 =- (76 +2116 | 
1d (26) 
=— = (16)=0. (22-32) 
From (22-32): 
£ 29 <9 (22-33) 
dt 


Integrating, we have r’6=constant. Denoting this constant by h, we have 
rO@=h (22-34) 
which is identical with (22-15),thus reaffirming the conservation of angular momentum. 


Defining a new variable u by 


ua ‘ or ee (22-35) 
r u 
and substituting in (22-34): 
6=hu’ (22-36) 
Now 
dt _ dr dé 
dt d6 dt 
-6# 94 (1) 
d@ dé \u 
Using (22-36) in the above equation : 
du 
| ft=—-h a0 (22-37) 
Also, 
di dr dé dr 


~ dt d@ dt ~ dé 


22-8 
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Using (22-36) and (22-37) in the above equation: 


2 
22d u 


f=—h'u 30? (22-38) 


Finally, using (22-35), (22-36), and (22-38) in (22-30): 


which is a linear differential equation equivalent to (22-30). 





+u=—— (22-39) 


22.3 CORRELATED VELOCITY 


This section is a detailed version of the method of determining the correlated velocity 
described in Reference (7),with the gaps and details filled in by the author of this report. 


There are several approaches to the solution of the problem. Although the method adopted 
here is not efficient in computation time, it was chosen because it is, from the author’s viewpoint, 
the least difficult and least prone to losing the insights into what is going on. The most efficient 
method may be somewhat difficult in providing insight for beginners. Those readers, who are 
interested in further study of the subject, are referred to an excellent book by Richard H. Battin 
entitled “ An Introduction to the Mathematics and Methods of Astrodynamics” (AIAA Education 
Series, 1987.), which must be counted among the great classics in scientific literature as the Editor- 
in-Chief of the series noted. 


Recalling that (Figure 22-1) 


h=r,xi=r,V,siny (22-40) 


and referring to the S term in (22-39), we may write 





BH. 
h? r,V,sin’y 
_ 1 
7 2 
Ty Vo oi 2 
r,sin’Y 
7 Si (22-41) 
Ar, sin’ y 


22-9 
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where 
r,V- 
14 6 6 «0 
L 
Substituting (22-41) into (22-39) : 
d’u(@ 1 
+u(§}=————__ 
do. (8) Ar,sin’y 


Referring to Figure 22-1, assuming 8=0 when t=0 andr(0) 
necessary for the solution of (22-43) are: 


and referring to Figure 22-2, 





(22-42) 


(22-43) 


=r, , the initial conditions 


(22-44) 





FIGURE 22-2 DEPICTION OF Cot y = 


d qd 1 
a9 29) dé r(@) 
__1 dr 

r? dé 
_ i dt 
x rdO 

22-10 


rd@ 


(ee _siieseaicee, Ewes Siprereres, ine, i». | 
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ee! 
ee coty 
Thus, since r=1rq when 0=6): 
du (6 1 cosy | 
= 22-4 
dO leo Tp Siny = 





The complete solution of (22-43) is (as shown in the appendix at the end of this section): 


1 1 (1—cos® | sin(y-9} 
u(6)=——__ = —_ | —————- + 22-46 
(8) r(@) Ty | Asin’y  —s- SY : 


The validy of (22-46) as the solution of (22-43) may be confirmed by differentiating u (6) W.R.T 


3 








du (6 _1 sinO cos(y-8 (22-47) 
dé Ty | Asin sin'y 
and 
2 
Pee aeedee Q ) (22-48) 
dé Ty | Asin“y siny 
and substituting (22-46) and (22-48) into (22-43). 
When @ reaches » (range angle), r becomes ry, the target vector. 
Evaluating (22-46) at this boundry condition, 
r ms in(y—9) | 
“o8 1 — re sin (Y b (22-49) 
Tr | Asin“y siny 
Solving (22-49) for 4 and recalling the definition of (22-42) : 
r,V- 1 
he c ° = (22-50) 
int sin (o —y )siny 
T 


Solving (50) for Vo: | 
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1 
1—cosd 2 


V.= E (22-51) 


. 4 : 2 
si ¥ +sin(o —y ) siny 


0 
.T 


There are many combinations of Vg and y (magnitude and direction of correlated velocity) 
which satisfy (22-51) for given values of ro, rr and . To determine the unique values of Vo and 
Y, we have to specify the time of (free fall) flight, which is shown next. : 


Referring to (22-16), integrating dt from 0 to t, (time of flight) and d@ from O tog (range 
angle), we have 


—f* 


(22-52) 


oh, 
e 
tl 
a oe 
ot, & 
a 
ww 
cD 
OQ. 
cD 





i i —_ -2 
—— [ae 8 } a6 (22-53) 


Asin’ y siny 


with A given by (22-50), and Vo given by (22-51). 


Equation (22-53) shows that for the given time of flight te and range angle $, the Vo and jy, 


which satisfy (22-53) becomes the correlated velocity V,. = Vo with the flight pass angle y, = at 
t = 0 or at deployment. 


The integrand of (22-53) may be written 


1-cos@ | sin is 0) | atic cos0 sinycos@ — cosy sin® F 
Asin’ y siny Asin’y Asin’y siny 











~2 
-| ~ +j1- ; cos8—corsin® 
Asin’ y Asin” y 
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=[a+bcos@+csin@]~7 (22-54) 


where 





Using (22-54) in (22-53) : 





d 
r 
pees | eee aims (22-55) 
9 (at bcos® + csin6 ) 


The integration may be performed, resulting in a closed form solution based on the formula table. 


After some algebra, the final algebraic equation for ts may be written : 


t, (1—cosd )coty + (1—A) sind ¥ 2siny a _ ((2/a)—1)? (22-56) 


VV. sin r 2 
_ (2-2) _ 
T 


3/2 6 
X x —] siny cot ($ |- cosy 
Referring to (22-56) and equivalently to (22-53), we note that tp becomes infinite at y =0, 
since siny is in the denominator. However, for the current LRBM applications, ‘y never becomes 
zero as explained below. 





Multiplying both sides of (22-49) by Asin”, we have 


Ir 
me Asin?y = 1—coso + Asinysin (y — 6 ) (22-57) 
T 


It follows from setting y=0 in (22-57) that we will have cosb = 1 or the range angle o = 
Q. This means that the current position vector is rt colinear with the target vector r, which never 
happens in the existing LRBM practices. Of course, this is intuitively obvious without 
mathematics. 


The objective of the following 1s to verify Equation (51) for the special case of a circular 
orbit (satellite injection into a circular orbit). 


22-13 
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For a special case of a circular orbit, |r, | =| 1,1! and y=90°. From the range angle 6 


introductory physics, we know, for a circular orbit of ro that 


w 





pm ae 
r Ty 
It follows: 
HL 
V,= 


Evaluating (22-51) with ry =r; and y= 90" and noting that 
sin (® —'y )=sind cos90 ° — cosd sin90 ° 
=—cosd 
we have 


_| p 1-cosd 
v.-|# 1—cosd 


=f 
ote 


Thus, verifying the validity of (22-51) for the special case of a circular orbit. 





Now, for a circular orbit of period Tp, V,T,=2nr,. That is 





For a range angle » and the corresponding time of flight t, in a circular orbit : 


t 


a8 
T, 21 
Using (22-62) in (22-63), 


22-14 


(22-58) 


(22-59) 


(22-60) 


(22-61) 


(22-62) 


(22-63) 
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ti (22-64) 








7 =] (22-65) 
Using y= 90°, A=1 in (22-53) and noting 


sin (y —@)=sin90 ° cos® —cos90 ° sin® 


=cos0 (22-66) 
we have 
> 
_ To -2 
rate || = a +cos@ eres | ay 
0 
aa ) (22-67) 
which is equal to (64). 


Thus verifying the validity of (53) for the special case of a circular orbit. 
22.4 TRAJECTORY: A SEGMENT OF ELLIPTIC ORBIT 


The general solution of (22-39) is given by 





u(@)= : = (1+ec0s0) | (22-68) 


mt 
@® 

= 
w 


arr esin® (22-69) 
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sl (6 =— ecos6 (22-70) 
dé h 


and substituting (22-68) and (22-70) into (22-39), and obtaining the right-hand side of (22-39). 


Rewriting (22-68) in a more conventional form: 


+(@)=—b_/_ (22-71) 
~ 1+ecos0 
h2 
In (22-71), “=a (1—e7) as shown in Appendix B. Thus, (22-71) becomes : 
a (1 - e?| 
o> 1+ecos0 ere) 


which is the equation for the conic sections (see for instance, “Calculus and Analytic Geometry” by 
Thomas, Addison Wesley) where e represents the eccentricity of the conic sections, and a is the 
semi-major axis in the case of ellipses. It is well known that the conic section becomes a circle if e 
= (0, an ellipse if 0<e<1, a parabola if e=1 and a hyperbola if e>1. 

Now, solving (46) for r(8): 


- 2 
T,Asin’y 


r(8)}= 1+[Asiny (siny —@)—cos6 | eee 


which is in the form of (22-72). 


From (22-72) and (22-73), we can see that the relationships between e and A may be determined. 
It turns out, after some algebra, that : 


A=1 corresponds to e =Q for circle ( with y =>). 
Q0<A<2 corresponds to 0<e<1 for ellipse 

XA =2 corresponds to e = 1 for parabola 

XN >2 corresponds toe > 1 for hyperbola 


In applications for LRBM, the combinations of ro, vo and y is such that A is always in the range 


of 0<A<2, which corresponds to 0 <e < 1, thus assuring the free falling ballistic trajectory is 
always a segment of an elliptic orbit. 
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22.5 COMPUTATION 


The key equation for the computation of the correlated velocity is the Equation (22-56). 


_ What we want is the magnitude of the velocity Vo and its direction y which satisfy (22-56) for a 


specified time of flight t,, with the given values of the initial release or deployment position Ip, the 


target position rr, and the range angle  betweenr, and ry. Since by the defintion of the dot 
product 


I, ° Lp=Tyr,cosd 


we can determine 6 from 





r 
-} *0 
= COS 
Yr 


oe 
= 
aq 15° 


Note that A in (22-56) is given by the definition : 


V=—— 22-42 
7 ( ) 


This implies that Vo appears in several places in (22-56) either explicity or implicity. 


For the determination of Vo and y by a numerical method, we start with an initial trial guess of j, 
say 45°, and get the value of Vo corresponding to this value of y by means of equation (22-51). 


Next, substitute the values of y and Vo obtained above to see if the right-hand side of (22-56) 
matches the LHS of (22-56), which is the specified value of the time-of-flight tp. Otherwise, we 


continue the process by adjusting the trial value of ‘y, which yields a new value of Vo, until the 
difference between the LHS and the RHS of (22-56) falls within the acceptable limit. The Vo and 


Y, which satisfy this condition is the correlated velocity we are looking for under the specified 
conditions. 
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APPENDIX A 
SOLUTION OF EQUATION (22-43) 


A-1/A-2 
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Repeating (43): 
d*u(0 gies 
de AV ,sin’y 
Thus, with = =D 
> d9 _ > 
: 1 
(D?+1)u=———— = K (A-1) 
Ar, sin” y 


Let U;, denote the homogeneous solution and U, particular solution of (A-1). Then for the 
homogeneous solution: 


(D?+1)=0 
(D+i)(D-i)=0 
D=i D=-i 

Thus, 

u, = C,cos0+ C, sin@ 
For the particular solution: 
Assume 
Up=A 
then 


u' =0 
P 


iI 
u 
P 


0 


Substituting into (A-1) : 
O+A=K 


so, the complete solution u=u, + u, iS 


u= C,cos6+C,sin6 +k (A-2) 


Repeating (22-44) 





u(0)=— at 8=0 
0 
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Substituting (A-2) evaluated at 0 = 0 into (A-3): 


+ =C +k 
Ir 


0 


or 


oat 


0 
Differentiating (A-2): 


a =-C, sin8 + C,cos0. 


Repeating (22-45): 


du 1 cosy 


a0 eo © ‘Siny 





Substituting (A-6) into (A-5) and evaluating at @=0: 


__1 cosy 


2 * 
T, siny 





Substituting (A-4) and (A-7) into (A-2) with K=———— 


cosy . 
Y in® + 





1 
u=|— -———_ |cos8- —- ——~s 
T Ar,sin’y YT siny 








(A-3) 
(A-4) 
(A-5) 
(A-6) 
(A-7) 

1 

Ar, sin ony 
Ar, sin oy 





A-4 


cxemsttomitatimie, els = eee ee ee eee eee Cee 
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u=-— 


1 |(1-—cos0 : siny cos8 — cosy sin® 
To {| Asiny siny 


— , Sin (y -8) | 


Ty | Asin’y siny 


which is identical with (22-46) in the main text. 


A-5 
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APPENDIX B 
EQUATION OF ELLIPSE 


B-1/B-2 
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FIGURE B-1. ELLIPSE 
The equation of an ellipse in rectangular coordinate is, as we all know, from analytic geometry: 


x y’ 
—+-—=1 (B-1) 
ab’ 
The equivalent form of an ellipse in polar coordinates with r(@) measured from one of the focal 
points (f) is given by : 


_a(i-e?} 
~ 1+ecos6 


r(@) 


(B-2) 


See, for instance, Thomas: Calculus and Analytic Geometry) where the eccentricity e is the ratio 


f ee : 
of distance (of) to (0a) i.e. e= - | The e has to be positive and less than 1 (i.e., 0 <e < 1) to 


have geometrical shape of ellipse. 


The semimajor axis (a) may be chosen arbitrarily as long as it is a positive real number. So, we 
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pick a number such that 
aw" | 
1-e’ — (B-3) 
where 
h= = rxf 
is the angular momentum per unit mass and 
L=GM, 


in which G is the universal gravitational constant, and Mg is the mass of the earth. 


Since 0 < e < 1 for an ellipse, e2 < 1. Thus, the RHS of equation (B-3) assures a > 0 because h? > 
Oandu>0. 


Substituting (B-3) into (B-2), we have 


h?/p 
r(8)= 1+ecos8 


which is identical with (22-71). 
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